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Abstract 



Quantum Chromodynamics with two light quark flavors is considered in the lat- 
tice regularization with improved Wilson fermions. In this formulation chiral 
symmetry is explicitly broken by cutoff effects linear in the lattice spacing a. As 
a consequence the isovector axial currents require improvement (in the Symanzik 
sense) as well as a finite renormalization if they are to satisfy the continuum 
Ward-Takahashi identities associated with the isovector chiral symmetries up to 
small lattice corrections of O(a^). 

In exploratory numerical simulations of the lattice theory algorithmic difficul- 
ties were encountered at coarse lattice spacings. There the hybrid Monte Carlo 
algorithm used suffers from a distorted Dirac spectrum in the form of unphys- 
ically small eigenvalues. This is shown to be a cutoff effect, which disappears 
rapidly as the lattice spacing is decreased. An alternative algorithm, the poly- 
nomial hybrid Monte Carlo algorithm, is found to perform significantly better in 
the presence of exceptionally small eigenvalues. 

Extending previously used methods both the improvement and the renor- 
malization of the axial current are implemented non-perturbatively in terms of 
correlation functions formulated in the framework of the Schrodinger functional. 
In both cases this is achieved by enforcing continuum Ward identities at finite 
lattice spacing. Together, this restores the isovector chiral symmetry to quadratic 
order in the lattice spacing. With little additional effort the normalization factor 
of the local vector current is also obtained. 

The methods developed and implemented here can easily be applied to other 
actions formulated in the Schrodinger functional framework. This includes im- 
proved gauge actions as well as theories with more than two dynamical quark 
flavors. 
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Zusammenfassung 



Quantenchromodynamik mit zwei leichten Quarks wird in der Gitterregular- 
isierung mit verbesserten Wilson Fermionen betrachtet. Die chirale Symmetrie in 
dieser Formulierung wird von Gitterartefakten, die linear im Gitterabstand a sind, 
explizit gebrochen. Daher erfordern die axialen Isospin Strome Verbesserung (im 
Symanzik Sinn), sowie eine endliche Renormierung, wenn sie die Ward-Takahashi 
Identitaten des Kontinuums bis auf kleine Gitterkorrekturen proportional zu 
erflillen sollen. 

Algorithmische Probleme bei groBen Gitterabstanden machen die numerischen 
Simulationen der Gittertheorie schwierig. Der Hybrid Monte Carlo Algorithmus 
leidet unter einem verformten Dirac Spektrum in Form unphysikalisch kleiner 
Eigenwerte. Es wird gezeigt, daB dies ein Gitterartefakt ist, welches schnell ver- 
schwindet, wenn der Gitterabstand verringert wird. Ein alternativer Algorith- 
mus, der polynomische Hybrid Monte Carlo Algorithmus, zeigt erheblich bessere 
Eigenschaften im Umgang mit den auBergewohnlich kleinen Eigenwerten. 

Durch Erweiterung und Verbesserung vorher verwendeter Methoden wird die 
nicht-perturbative Verbesserung und Renormierung des Axialstroms durch Kor- 
relationsfunktionen im Schrodinger Funktional implementiert. In beiden Fallen 
wird dies erzielt, indem man Ward Identitaten des Kontinuums bei endlichem 
Gitterabstand erzwingt. Zusammen stellt dies die chirale Symmetrie bis zur 
quadratischen Ordnung im Gitterabstand wieder her. Mit wenig zusatzlichem 
Aufwand wird auch der Normierungsfaktor des lokalen Vektorstroms berechnet. 

Die Methoden, die hier entwickelt und implementiert wurden, konnen leicht 
auch flir andere Wirkungen verwendet werden, die im Schrodinger Funktional 
formuliert werden konnen. Dies umfaBt verbesserte Eichwirkungen sowie Theo- 
rien mit mehr als zwei dynamischen Quarks. 
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Chapter 1 
Introduction 



Compared to the electroweak sector of the standard model of particle physics, 
quantum chromodynamics (QCD) with its few parameters and extensive symme- 
tries seems to be a rather simple theory. Still, it is expected to describe the whole 
spectrum of strong interaction phenomena from high-energy hadron collisions to 
the decays of heavy-quark bound states and of course the hadron masses them- 
selves. Since its birth in the 1960s QCD has been confronted with experiment in 
innumerable cases and is now a well-established part of the standard model. 

At high energies the relevant degrees of freedom (quarks and gluons) are found 
to be only weakly coupled and the non-interacting theory can serve as a starting 
point for a perturbative treatment. The breaking of Bjorken scaling [1] in deep 
inelastic lepton-hadron scattering, which is the original and still one of the most 
powerful quantitative tests of (perturbative) QCD, is associated with these energy 
scales. 

At low energies the strong interactions show a very different behavior. The 
coupling becomes strong such that a description in terms of weakly interacting 
quarks and gluons is no longer appropriate. Instead, the relevant degrees of 
freedom seem to be the light mesons (pions). Also here precise experimental 
data are available, in particular for the masses and other properties of the light 
hadrons. If one wants to establish that the same Lagrangian describes both 
regimes, vast energy differences need to be bridged [2]. 

As an additional complication, the transition in the effective degrees of free- 
dom makes it much more difficult to work out the QCD predictions at low en- 
ergies, because perturbation theory, the "standard tool" of particle physics, fails 
here. While effective theories are useful at this point, they do not constitute a 
first-principle method. Also, some non-perturbative predictions can be made on 
the basis of symmetry considerations or QCD sum rules, but these can not be ap- 
plied to all problems one might be interested in and often also require additional 
assumptions. 

With the lattice regularization Wilson [3] proposed a radically different ap- 
proach. Lattice QCD allows for first-principle predictions without any additional 
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assumptions and can be applied to a variety of problems, from the hadron spec- 
trum to meson decay constants and the running of the coupling or even the 
topological structure of the QCD vacuum. 

In this formulation one discretizes the fields and the action using a Euclidean 
space-time lattice. Thus, before any computation is made, the high frequencies 
are removed from the theory, rendering the latter ultraviolet finite. The infrared 
divergencies related to zero modes can also be be cured by either considering 
a non-zero quark mass or a finite volume with fixed boundary conditions in at 
least one direction (see Chapter EI). At this point the theory is mathematically 
well-defined without reference to perturbation theory. 

Moreover, if the lattice is small enough, a numerical evaluation on the com- 
puter becomes feasible. In this way the lattice serves two purposes: it regularizes 
the theory with a momentum cutoff proportional to the inverse lattice spacing 
and at the same time it is a tool to evaluate observables non-perturbatively. 
While today the main aspect of lattice QCD is clearly in numerically obtaining 
phenomenological predictions from Monte Carlo simulations, one should not for- 
get that in many cases the lattice is the only way to define a quantum field theory 
beyond perturbation theory. 

As with any other regularization, the regulator has to be removed before 
results can be compared to experiment. On the lattice this means that the 
continuum limit has to be taken by making the lattice spacing a smaller and 
smaller. Here one encounters the usual ultraviolet divergencies, which require 
renormalization of the bare parameters and operators. 

In numerical simulations the accessible lattice spacings are strongly con- 
strained by the computer resources and the scaling behavior of the available 
algorithms. In fact, for a long time after the invention of lattice QCD, simula- 
tions were only possible in the quenched approximation. Here the observables are 
evaluated on a gauge background generated with the gauge action only, which in 
terms of Feynman diagrams amounts to removing all virtual quark lines. While 
this makes the simulations significantly easier, it does not represent a controlled 
approximation and thus (in principle) calls into question all physical predictions 
obtained in this way. Still, impressive results were obtained from quenched sim- 
ulations and in particular the hadron spectrum is in rather good agreement with 
experiment [4]. 

With new computer generations and algorithmic improvements unquenched 
(or " dynamical" ) simulations are now possible at reasonable lattice spacings and 
volumes. However, these algorithms typically slow down proportional to a'' (see 
e.g. [5]) or even worse, such that a factor two in lattice spacing can change 
the computational effort by more than two orders of magnitude. As a result a 
is usually varied only in a very limited range, say from 0.1 fm to 0.05 fm and 
hence close attention should be paid to how the observables of interest approach 
the continuum. Naturally, this is strongly influenced by the details of how the 
continuum theory was discretized. 
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Today there exists a large variety of discretizations of lattice QCD in addition 
to Wilson's original formulation, on which this work is based. These differ in both 
the discretization of the gauge and - more importantly - the fermionic action. 
A brief overview with emphasis on their different properties concerning chirality 
will be given in Chapter |3J 

A systematic description of the approach to the continuum limit of a lattice 
theory was found by Symanzik [6-8]. The symmetries of the lattice action define 
a set of operators, which, if inserted in correlation functions of the continuum 
theory, asymptotically reproduce the discretization errors of the lattice theory. 
In other words, the lattice theory is described in terms of an effective low-energy 
(•C a~^) continuum theory, whose Lagrangian contains (4 + fc) -dimensional in- 
teraction terms with couplings proportional to a^. On this basis one exploits the 
freedom to change the lattice discretization by irrelevant operators to remove the 
cutoff effects, i.e. the lattice artifacts, order by order in the lattice spacing. 

This procedure is known as the Symanzik improvement programme and con- 
cerns both the action and the composite fields, i.e. for QCD in particular the 
quark bilinears. In this way the continuum limit can be "accelerated" and a 
smaller range of lattice spacings might be sufficient for a reliable extrapolation. 
The application of the improvement programme to Wilson fermions is a central 
part of this work and a first (general) discussion will be given in Chapter El 

Renormalization is commonly discussed in the framework of perturbation 
theory, where divergencies in Feynman diagrams are removed e.g. through a 
subtraction prescription based on dimensional regularization. Due to the non- 
perturbative nature of the low-energy sector of QCD this approach is no longer 
sufficient in this case. As will be detailed in Chapter |H1 a perturbative calculation 
of renormalization factors does not result in a controlled estimate of systematic 
errors of renormalized correlation functions at low energies. 

When approaching the continuum, the bare parameters in the action have to 
be tuned such that a set of renormalized quantities is kept fixed. These can be 
mass ratios or otherwise defined renormalized couplings and their specific choice 
defines the renormalization scheme. In this context as well as for the additional 
renormalizations required for the composite fields we employ a finite volume 
scheme based on the Schrodinger functional first discussed in [9,10]. Chapter El 
will provide the reader with the necessary background. 

Despite its rich structure, some of the phenomenologically most relevant ques- 
tions are not concerned with effects within QCD alone. In particular, one often 
considers QCD observables, which can be interpreted as matrix elements of the 
effective weak Hamiltonian between QCD bound states. A prominent example is 
the pion decay constant F^, defined through the matrix element 



of the axial current between the vacuum and a pion state with momentum p. 
Regarded as an insertion of the effective weak Hamiltonian, it parameterizes the 



{0\A^{0)\7i) = ip^F^ 
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amplitude for a pion to decay into a lepton anti-neutrino pair through a virtual 
W boson. Here it is most evident how the renormalization and improvement of 
the axial current directly affect the physical result of a lattice estimate of F^. 

For quenched lattice QCD with Wilson fermions the renormalization fac- 
tor [11] and improvement constant [12] of the axial current are known non- 
perturbatively. Both are obtained by enforcing the chiral symmetry of the un- 
derlying continuum theory at finite lattice spacing. The non-perturbative im- 
plementation for the case of lattice QCD with two degenerate flavors of Wilson 
quarks is the central part and the resulting renormalization factor and improve- 
ment coefficient are the main result of this work. 

Another important application of the renormalized improved axial current is 
in the context of a programme to calculate the fundamental parameters of QCD 
from hadronic input parameters. In fact, the work presented here is an integral 
part of an effort to calculate renormalized quark masses from first principles. The 
former also constitute an essential ingredient in tests of the standard model. 

With the same method as in the quenched case [13], the starting point for such 
a calculation is the current quark mass derived from the non-conservation of the 
axial current (PCAC relation). The renormalization of the pseudo-scalar density, 
which is also required in this context, has already been performed [14, 15] such 
that together with the results from this work all necessary tools are available. 

Starting from a review of the most important properties of QCD (Chap- 
ter |2|), which includes a derivation of the Ward-Takahashi identities associated 
with the chiral symmetry, we will turn to the peculiarities of the lattice formula- 
tion (Chapterini) and discuss the continuum limit, renormalization and Symanzik 
improvement. This is followed by a chapter concerned with the status of chiral 
symmetry in different lattice discretizations (Chapter|3]) with particular emphasis 
on Wilson fermions. The introductory part concludes with a presentation of the 
Schrodinger functional (SF) as a renormalization scheme and our chosen method 
to compute renormalization factors and improvement constants (Chapter Ej). The 
necessary correlation functions and notation to be used are also introduced there. 

The first chapter of the main part (Chapter ^ discusses Monte Carlo algo- 
rithms, data analysis and other technical aspects of our simulations. Algorithmic 
issues, which we faced in the numerical evaluation of the axial current normaliza- 
tion condition, are also reported. These problems are traced back to a distortion 
of the spectrum of the Wilson-Dirac operator, which in turn can be interpreted 
as a cutoff effect. 

The remaining chapters are devoted to the non-perturbative axial current 
improvement (Chapter Ej) and renormalization (Chapter IH}. The former is im- 
plemented by requiring a current quark mass derived from a Ward identity to be 
independent of the external states, which are varied using projection techniques. 
The integrated axial Ward identity with operator insertions is used to formulate 
a normalization condition for the axial current on the lattice. Numerical results 
of our implementation are presented and summarized in interpolating formulae 
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for future use. 

Finally, the last chapter summarizes all results and discusses their possible 
application. In abbreviated form the findings of the present work have been 
pubhshed in [16-20]. 



Chapter 2 
Continuum QCD 



2.1 History and properties 

Already in the 1960s it was conjectured that the observed large number of strongly 
interacting particles, the hadrons, are composite objects made from supposedly 
fundamental building blocks called partons. In the spirit of Rutherford's scatter- 
ing experiments, which revealed the structure of the atom, hadrons were probed 
with beams of highly energetic leptons in order to learn about the underlying 
dynamics that govern the formation of partons into hadrons. 

The most interesting experimental results came from the kinematic region of 
deep inelastic scattering, where both the momentum transfer and the energy 
transfer u from the leptons are very large with the ratio q'^/iy fixed. As proposed 
by Bjorken in 1969 [1], in this region the structure functions, which parameterize 
the momentum distribution within the hadron, were found to depend only on the 
ratio q'^/jy (Bjorken scaling). 

The easiest way to understand this behavior is to assume that the leptons 
scatter off almost-free pointlike particles, the constituents of the hadrons. To 
accommodate Bjorken scaling, the theory describing the dynamics of the partons 
should therefore have the feature that the interaction becomes weak at high 
energies (or small distances). The only known quantum field theories^ with this 
property, which is now called asymptotic freedom, are the non-Abelian gauge 
theories introduced by Yang and Mills [21]. 

The working hypothesis is now that the interaction of the fundamental degrees 
of freedom - at this point renamed quarks by Gell-Mann - is described by a non- 
Abelian gauge theory. 

At the same time experimental observations required the quarks to have an 
additional unobserved quantum number (" color" ) in order to avoid conflict with 
the Pauli exclusion principle. Most of the difficulties could be resolved by iden- 
tifying the symmetry corresponding to the new quantum numbers with the non- 

"'^More precisely: the only rcnormalizablc quantum field theories in four dimension. 
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Abelian gauge symmetry. Further experimental input and tlieoretical arguments 
uniquely fixed the gauge group to be SU{3). This theory was named quantum 
chromodynamics, QCD, and its gauge quanta are called gluons. 

At this point one should note that it is the property of being asymptotically 
free that ensures the applicability of perturbation theory at high energies. One of 
the successes of perturbative QCD is the correct prediction of the (logarithmic) 
corrections to Bjorken scaling, which is now known to be strictly valid only in 
the limit of infinite energy. 

In contrast to Abelian gauge theories like quantum electrodynamics, QED, 
the field quanta of non-Abelian theories also carry charge and thus interact 
among themselves. Indeed, perturbation theory shows that it is precisely this 
self-interaction, which is responsible for asymptotic freedom. The situation is 
opposite to the one found in QED, where perturbation theory predicts the effec- 
tive charge of the electron to increase with energy. 

We can now write down the Lagrangian density for the gluonic (Cg) and 
fermionic (£/) part of QCD with A^f flavors of quarks. 

with the fleld strength F^^ = d^A^ - d^A^ + g^lA^, AJ\ . 

Here go is the bare gauge coupling and A^ the Lie-algebra valued anti-Hermitian 
SU (3) gauge fleld. For the fermionic part we have 

I^S = (2.2) 

/=i 

with the gauge covariant derivative Ip = '~ff^{d^+goA^) and the bare quark masses 
rrif. The two color indices of as well as the color and Dirac indices of ip and 
ip have been suppressed. 

To account for the observed particle spectrum, hadron states and in general 
all physical observables are postulated to be color singlets. In this way an un- 
observed proliferation of states due to the color symmetry is prevented. This 
non-perturbative phenomenon (known as confinement of color) is ascribed to an 
increase of the force between color sources at long distances (~ Ifm). Since, as 
a consequence, even at high interaction energies the initial and flnal states are 
subject to conflnement, also perturbation theory is affected and hadronic matrix 
elements cannot be obtained using a purely perturbative treatment. This prob- 
lem is usually addressed by factorizing them into a "hard" (perturbative) and a 
"soft" (non-perturbative) part using e.g. the operator product expansion (OPE). 
The soft part is then parameterized by effective couplings. 

Lattice simulations of pure non-Abelian gauge theories show that for large 
distances the energy of two static color sources grows linearly with the so-called 
string tension a ^ IGeV/ fm (see e.g. [22]). Hence there is strong numerical 
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evidence in support of tlie confinement hypothesis. However, an analytical first- 
principle explanation from the underlying dynamics is still missing although many 
attempts have been made, e.g. in terms of condensation of topological excitations 
[23]. 

Finally, the most important ingredient to understanding the strong interac- 
tions at low energies is chiral symmetry. For degenerate quarks the QCD La- 
grangian ()2.2|1 is flavor-blind and hence invariant under (global, i.e. space-time 
independent) unitary transformations of the A'^f-component fermion fields. In 
addition, for massless quarks there is also no coupling between the left- and 
right-handed field components 

^L = |(1-75)V^ and = |(1+75)V^. (2.3) 
Therefore separate U{N{) transformations of these fields according to 



^ ^ Ur^r + Ul^l (2.4) 



are also a symmetry of the massless QCD Lagrangian. These are usually written 
in terms of the vector and axial vector transformations 

Ur = Ul = U ^Uifj = e'^^^''^lj (vector) 

Ur = uI = U ^^U^il) = e*^"^''^^^ (axial) . ^'^'^^ 

At the classical level this U{N{)y x U{Ni)a symmetry gives rise to two conserved 
currents (the vector and the axial current) for each generator /° of U{Ni), 

where the 1°' act on the flavor indices of the fermion fields. To summarize, the 
concrete requirements for the U{N{)y x [/(iVf)A flavor chiral symmetry are that 
the Dirac operator is diagonal in flavor space and that it anti-commutes with 75. 

Since the masses of the up and down quarks are much below typical hadronic 
energy scales, massless QCD with iVf = 2 seems to be a reasonable approximation 
of nature. In this case, the vector symmetry SU{2)y x U{l)y implies isospin 
symmetry and quark number conservation, which are indeed experimentally con- 
firmed to a high precision. 

The fact that the corresponding axial symmetries are not found in the strong 
interactions is believed to be due to a condensation of quark-ant iquark pairs, 
such that these symmetries are spontaneously broken by the QCD vacuum. In 
this picture the isospin triplet of light pseudo-scalar mesons, the pions, become 
the quasi-Goldstone bosons of the spontaneously broken axial symmetry.^ Their 



direct proof of the Goldstone theorem for this case will be given in Section UTTl 



2.2 Euclidean path integral 
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non-vanishing mass is explained by the fact that for small but non-zero quark 
mass the SU{2)x is only an approximate symmetry and the pions are hence not 
true Goldstone bosons. 

The identification of the pions as the Goldstone bosons of spontaneous chiral 
symmetry breaking leads to various relations between matrix elements, whose 
experimental verification helped in establishing this picture. A systematic study 
of the low-energy limit of QCD can be performed in the framework of chiral 
perturbation theory [24]. This is a low-energy effective theory of the strong in- 
teractions, where pions are introduced as fundamental degrees of freedom. The 
quark masses are treated as a perturbation to the chirally invariant Lagrangian. 

In a somewhat worse approximation of nature this picture can be extended 
to three flavors if the strange quark mass is also neglected. The Ward-Takahashi 
identities resulting from the SU (2)v x SU (2)a symmetry of massless two-flavor 
QCD are derived in Section ESI 

The conservation of the singlet axial current corresponding to the remaining 
U{1)a symmetry is violated by quantum corrections, which are related to the 
topological structure of the QCD vacuum. This is known as the axial (or chiral) 
anomaly and is crucial in explaining the absence of an isosinglet pseudo-scalar 
meson with mass comparable to that of the pions. In the path integral formu- 
lation (see below) the anomaly can be understood from the fact that the U{1)a 
transformation does not leave the integral measure invariant [25]. This is particu- 
larly transparent [26] in a lattice regularization using Ginsparg-Wilson fermions, 
which are invariant under a chiral symmetry, differing from ()2.5|) only at finite 
lattice spacing. 



2.2 Euclidean path integral 

While scalar field theories can be quantized efficiently in the operator language, 
for gauge theories it is more convenient to employ the functional integral (or path 
integral) formalism. In this approach the fundamental quantity is the generating 
functional 

Z[J] = I expli [ d^x{C + J- fields) I . (2.7) 

ificlds I i J 

Correlation functions are generated by taking functional derivatives with respect 
to the sources J(x) and perturbation theory is set up by expanding the (gauge 
fixed) action S = J dx^£ around a saddle point. Feynman rules can essentially 
be read directly from the Lagrangian. 

However, the functional integral is also the starting point for the so far most 
successful non-perturbative treatment of QCD, the lattice regularization. There 
the partition function Z = Z[0] is evaluated directly using Monte Carlo meth- 
ods. This approach requires a statistical interpretation of the path integral ()2.7|) . 
where the phase factor exp{iS) becomes a Boltzmann weight. To achieve this 
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the theory is transcribed from Minkowski to Euchdean space by introducing an 
imaginary time coordinate ("Wick rotation"). This formulation also reveals the 
deep connection between quantum field theories and statistical systems. 

With a purely imaginary time coordinate the space-time metric becomes Eu- 
chdean. Under certain condition the Minkowskian correlation functions (the so- 
called Wightman functions) can be analytically continued into this region to 
obtain Euclidean correlation functions (Schwinger functions). From the more ac- 
cessible Schwinger functions the physical Wightman functions in Minkowski space 
and thus the quantum field theory generated by ()2.7p can be reconstructed. 

In the path integral formulation the fermionic anti-commutation properties 
are accommodated through a representation of the fields ip and ip as Grassmann 
variables. The integral over the fields is therefore a Grassmann integral for the 
fermionic part. On the lattice the gauge fields are represented by compact link 
variables (see Chapter E)) such that the path integral is then defined in terms of 
the Haar measure on the gauge group. 

The QCD fermion action for one quark fiavor in its Euclidean form is given 

by 

S = jd^x^{i^D^+mo)^, , (2.8) 
where = d^+g^A^ . (2.9) 

Before the properties of the Euclidean lattice formulation as a non-perturbative 
regulator are discussed in Chapter El where for simplicity the case of a pure gauge 
theory is considered, we will derive the Ward identities associated with the local 
versions of the transformations ()2.5j) . 



2.3 Current algebra and continuum 
Ward identities 

For the two-fiavor case (A" = r"/2 with the Pauli matrices t"") the infinitesimal 
chiral transformations are written as 



vector <^ - - . 1 (2.10) 



We define the vector and axial variations of (composite) fields O through 
O -^0 + i60, with 



SyO = uj''6^0 and 6^0 = u^S^O , (2.12) 



2.3 Current algebra and continuum Ward identities 
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where (in case of a local transformation) uj"' is evaluated at the space-time point 
where O lives. For the quark fields we obtain 



2 



(2.13) 
(2.14) 



These can be used to calculate the variation of arbitrary expressions O built 
from the basic fields by treating 5" and 6^ as first order differential operators. In 
particular, for the variations of the isospin vector and axial vector currents, 



V^ix) = ^|Jixh,lT'^^^J{x) 



v4;(x) = V^(x)7^75|rXx 



(2.15) 



one obtains for example 



since the generators of SU (2) satisfy the algebra 



2 2 



2 2 



■abcTf_ 
Lt 2 



Calculating all variations using {7/^,75} = it is easily verified that the currents 
form a closed algebra 



-ie'^*'=y;(x). 



5eA^(x) = -ze'^''^A^(x), 5,M^(x) = -ie^'^V;;{x) . 



(2.16) 
(2.17) 



The Ward identities associated with the chiral symmetry of the action are now 
derived by performing local infinitesimal transformations of the quark and anti- 
quark fields in the Euclidean functional integral. 



2.3.1 Variation of the action 

The starting point is the Euclidean fermion action ()2.8j) for two quark flavors, 
where now niQ is a diagonal 2x2 matrix in isospin space. We perform local in- 
finitesimal variations of the fermionic fields and study their effect on the classical 
QCD action. These variations are parameterized by a;"(x) with support in a 
space-time region TZ. The isospin vector variation of the fermionic action is 



d^x 



n 



d^x 



n 



- t^>ir''(7^9^+mo)^ + ^(7^9^+mo)cj''|r''^ 
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where the terms containing the gauge fields canceled. Performing a partial 
integration of the second expression, which does not give boundary terms since 
the variation vanishes outside the region TZ, yields 



n 



iplr^l'j^d^ + mo]ip - 'Ji,{d^tp) - mo 



n 



d'^XLu'^{x) — dfj,V^{x) — ip{x) [|r", mo] 



(2.18) 



In the case of degenerate quarks mo is proportional to the unit matrix and hence 
the commutator with vanishes. 

In the same way the isospin axial vector variation of the action can be com- 
puted using the anti-commutator {7^,75} = 0. Since in ()2.14p the variations of 
ip and appear with the same sign, the mass now results in an anti-commutator 
term. Again performing a partial integration of the second term results in 



Em 



d^x 



n 



n 



d!^xuo°-{x) — 9^A"(x) + iI){x)'^^{\t°' ,mQ}%l){x) 



For two degenerate quarks this expression becomes 



5a5 



d xijj"'{x) 



n 



- d^Al{x) + 2moP(a;) 



(2.19) 



with the pseudo-scalar density 



P'^(x)=^(x)75|rXx). 



(2.20) 



2.3.2 Ward identities 

Through a formal manipulation of the functional integral the expectation value 
of the variation of an operator can be related to that of the action. Later we will 
discuss how such relations are realized in the regularized quantum theory. If the 
fermion integration measure is denoted by I^^, we have for a linear transformation 
of the Grassmann fields ^ 

^' = Aijj ijj' = ipA 
=^ [ViiVil)] = detAdetA[V^lj'Vip'] = J-^[ViJj'V^l;']. (2.21) 

^The rules of Grassmann integration imply that the Jacobian is the inverse of what one 
would expect from bosonic integrals, hence the appearance of J^^ in H2.21|l . 



2.3 Current algebra and continuum Ward identities 
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Using det(l+u;X) = l+tuTr X+0(c<j^), one obtains for the infinitesimal versions 
of the isovector transformations ()2.1()j) and ()2.1H1 

vector J = l + 0(u;2) (2.22) 
axial J = l + icj"Tr(r"75) + 0(cj2) . (2.23) 

Since the Pauli matrices are traceless, the Jacobian is unchanged for both isovec- 
tor transformations and one can perform the infinitesimal variable transformation 
in the path integral according to 

Z= I e-^ (C) = Z-^ I Oe~^ = I {0 + 50)e-^{l-5S) 

J fields J fields J fields 

= (C) + - (pbS) 
(SO) = {06S) . (2.24) 

We can now derive the integrated Ward identities associated with the flavor chiral 
symmetry of the theory with two degenerate quarks. To this end we insert the 
vector ()2.18j) and axial ()2.19j) variations into ()2.24j) . Choosing there an operator 
C^ext (for external) with support only outside the region TZ, where we perform the 
transformations, the left-hand side of eq. ()2.24|) vanishes and we have 



= Jji^xu''{x)^-d,y^{x)0,^t) (2.25) 
and = y"dScj''(x)^[-9^A;j(x) + 2moP(x)]Ccxt) • (2.26) 

These hold for any variation uj"'{x) and one can thus conclude that the expectation 
value multiplying it vanishes. This gives the vector current conservation in the 
form 

{d,V^{x)0,,,) = , (2.27) 

and for the isovector axial current one has the PCAC (partial conservation of the 
axial currents) relation 

{d,A';{x)0,,,) = 2m,{P\x)0,^,) . (2.28) 

Both are valid for operators Cext not located at the point x. Integrated Ward 
identities with non-trivial operator insertions will be derived in Section 18.11 to 
be used in the normalization conditions for the vector and axial currents in the 
lattice regularization with Wilson fermions. 



2.3.3 Anomalous symmetries 



Since the quantized theory requires a regularization, it is not immediately obvious 
how e.g. the above relations have to be interpreted. If a symmetry is preserved 
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in the quantum tlieory, this implies that it is either unbroken by the regulator or 
is recovered in the hmit where the latter is removed. Perturbation theory shows 
that this is the case for the isovector symmetries discussed above. In Chapter |S1 
we will see that in a specific lattice formulation of the theory the above relations 
can be made to hold up to (small) lattice corrections. This can be interpreted 
as showing that even beyond perturbation theory the isovector symmetries are 
preserved in the quantized theory. 

A transformation is called anomalous if it is a symmetry of the classical theory 
(i.e. the unregulated action), but not of the quantum theory. Depending on the 
chosen regularization, an anomaly can arise in different ways. One possibility is 
that the integration measure of the functional integral is not invariant under the 
transformation. If such a transformation is parameterized by an infinitesimal a;, 
we can write the Jacobian as J = 1 + iuj5J and instead of ()2.24|1 one then has 

{50) = {06S) + {06J) . (2.29) 

It can also happen that a symmetry is broken by the regulator in such a way 
that it is not recovered in the renormalized theory. To illustrate these points we 
consider the isosinglet transformations 

vector 6yxlj = ^lj , Sy^j = -^fj , 
axial 6j^ip = 75?/^ , (5a = ^^75 . 

For those eq. ()2.22|) is still true and in fact the U{l)y (quark number conservation) 
is preserved in the quantized theory. For the axial transformation, instead of 
fl2.2H|l we obtain for the Jacobian 6 J = 2Tr75. 

The axial anomaly mentioned in the introduction has the non-conservation 
of the isosinglet axial current = ip'jfj^'j^ip even at vanishing quark mass as 
a physical consequence. Therefore any admissible regularization of QCD needs 
to be able to reproduce it correctly. This can happen through a non-trivial 
redefinition of 75, i.e. for (j2.3()j) to be a symmetry of the regularized action one 
has to replace 75 with some 75. This is the case for dimensional regularization or 
Ginsparg-Wilson fermions (see Section l^!T|) . where the anomaly comes from 6 J 
in (j2.29|) according to (m = 0) 

(9^A^(x)axt) = 2(axtTr75) . (2.31) 

Alternatively, the anomaly can arise from an explicit chiral symmetry breaking 
by the regulator, such that the isovector axial symmetry is recovered in the 
renormalized theory, but the U{1)a is not. This the case in a lattice regularization 
with Wilson fermions, where the anomaly can be understood as coming from the 
chiral variation of the Wilson term (see Section 1^?^ . 



Chapter 3 

Gauge theories on the lattice 



Consider a field theory with a matter field ip^{x), where x is the space-time 
coordinate and i refers to some internal degree of freedom. The requirement 
that one should be able to choose the basis to describe this internal degree of 
freedom independently on every space-time point results in the principle of local 
gauge invariance. To achieve this, the matter fields need to be coupled to a 
gauge field A^^{x), which takes care of the basis transformation between different 
(but infinitesimally close) space-time points ("parallel transporter"). This is 
done via the "minimal coupling prescription", which replaces the derivative 
with the covariant derivative = {d^ + goA^) such that D^ip has the same 
properties under local basis transformations {gauge transformations) as ip itself. 
For the case of QCD the internal degree of freedom of the quark fields ip is that 
of "color" and i is the index of the fundamental (3) representation of SU{3) [ip 
is in the 3 representation). The anti-Hermitian gauge fields are conventionally 
parameterized as 

4 = E^Am(^")^.' (3-1) 

a 

where are the eight generators of SU{3) (Gell-Mann matrices). The trace of 
the square of the local field strength ()2.1|) is the only term that can be made from 
the gauge fields that is gauge-invariant, has mass-dimension four or less and is 
compatible with time-reversal, charge-conjugation and parity. The geometric 
interpretation of the field strength tensor F^iy is the local curvature in color space 
as obtained by connecting color bases around an infinitesimal square in the fi-u- 
plane. 

In the discrete theory the matter fields are defined on the sites of the lattice 
and to gauge connect those we always need parallel transporters for finite dis- 
tances. It is therefore more natural to describe the gauge degrees of freedom in 
terms of link variables Uf^{x) G SU{3). These are the path-ordered exponentials 
of the integral of ()3.1|) from x along the link of length a to the next lattice site in 
the n direction. The variable U^{x) is thus associated with this link and by defi- 
nition the link in the opposite direction is its inverse U-^{x+fi) = Ujj^{x), where 
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/i is a vector of length a in the // direction. Under a local gauge transformation 
L(a;) G SU{3) the links and quark fields transform as 

U^{x) L(x)f/^(x)L(x+/i)^ 

^ L{x)i/j{x) (3.2) 

^/^(x) ^ '?/'(x)L(x)''" , 

such that 

V^4j{x) = ^[U^{x)4j{x+fi)-4j{x)] 

V^x) = '-[^ix)-Ulix-f,Mx-f,)] (3.3) 

V^^(x) = l(V^ + V;,)^(x) 

define covariant derivatives on the lattice. Normal lattice derivatives d^, d* and 
dfj, are defined by setting f/^ = 1 in eqs. ()3.3|1 . 

On the lattice the field strength is measured by tracing the gauge links around 
closed loops and hence the simplest admissible lattice gauge action is Wilson's 
plaquette action 

SAU] = \Y.^r {1 - U,} , (3.4) 

where p runs over all oriented elementary plaquettes (1x1 loops) and Up is the 
product of the gauge links around a plaquette 

Up{x) = U^;^,, = U^{x)-U,{x+fi)-Ul{x + v)-U,{x)^ , (3.5) 
Up{x) A k{x)Up{x)k{x)^ . (3.6) 

For gauge group SU (N) the bare gauge coupling is commonly expressed through 
j3 = 2N / . A small a expansion of the gauge links in terms of the gauge fields 
gives 



U.;,,u = 1 + a^F^, + -Fl, + ... (3.7) 

such that the sum over oriented plaquettes in ()3.4p reduces to the Yang-Mills 
action ()2.1|) . In fact this is true for any lattice action of the form ()3.4|) composed 
of closed loops of gauge links. Improved gauge actions (see Section 13. 3j) can thus 
be defined by adding more extended loops with suitably chosen coefficients [27,28] 
to 



3.1 Monte Carlo integration 

To obtain the expectation value of a gauge-invariant observable F\U] in the 
quantum field theory described by ()3.4|1 we need to evaluate the integral 



(F) = ^ JvU F[f/]e-^[^l with Z = jvU e"^!^! , 



(3.^ 



3.1 Monte Carlo integration 
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where T>U denotes the product of the SU (3) Haar measures for all links on the 
lattice. For a 4-dimensional lattice of linear extension L there are 4-(L/a)^ links, 
which in typical simulations is of order 10^ or larger. The only way to evaluate 
such high-dimensional integrals is using Monte Carlo techniques. 

To this end one generates gauge link configurations f/j that are distributed 
with a probability Pj oc P[/je~'^t^'l (importance sampling). The integral ()3.8|) is 
then statistically approximated by 

1 ^ 1 
(P) = 7^Em]+0{^). (3.9) 

i=l ^ 

Unless a global heatbath algorithm is available for the system one wants to sim- 
ulate, configurations are produced by an algorithm that updates a given config- 
uration to obtain the next one (Markov chain). Configurations generated in this 
way are not statistically independent and therefore autocorrelation effects need 
to be taken into account in the data analysis. 

For the case of pure gauge actions (or equivalently in quenched QCD) the local 
heatbath [29,30] is an efficient algorithm to generate correctly distributed ensem- 
bles of gauge configurations and lattice sizes up to L/a = O(IO^) are possible 
today. When the correlation lengths in the system become large, autocorrelation 
along the Markov chain increases and the algorithm rapidly becomes inefficient 
(critical slowing down). This can be countered by adding (micro canonical) over- 
relaxation updates [31] to the heatbath algorithm. 

If we want to consider full QCD, the fermionic part needs to be treated an- 
alytically because the Grassmann-valued quark fields cannot be handled in a 
computer simulation. Since the fermion action is bilinear in the quark fields, 
the Grassmann integral can be performed and for the partition function one gets 
{M = D + m) 

Z = JvUV^pVije-^o^^^-^^^^^^^ (3.10) 

oc JvU detM[f/]e-^«[^l . (3.11) 

In the path integral for an observable 0{U, ip, ijj) this integration of the fermionic 
fields results in a new observable 0{U), which depends only on the gauge field 
(Wick contraction). Instead of the fermion fields O contains matrix elements of 
the fermion propagator S = and as a consequence it can be a very compli- 
cated and non-local function of U. For the "pion-pion" correlator ('075'0)x('075'0)2/ 
one obtains e.g. 

((^75^)x(^75V^)y) = 

^ JvU (Tr{5..[f/]75}Tr{5,,[f/]75} 

-Tr {S,y[U]-f,SyMl5}) detM[f/]e-^«[^l , (3.12) 
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where the trace is taken over Dirac indices. Algorithms for simulating fuU QCD, 
i.e. the generation of a gauge ensemble according to VU detM[U]e~^''^'^\ will be 
discussed in detail in Chapter IHl 

3.2 Continuum limit and renormalization 

As already mentioned, the regulator has to be removed before physical predictions 
can be made. In the statistical analogy the mass (in lattice units) a-m of the 
lowest-lying state (e.g. the 0"*"+ glueball for a pure gauge theory) corresponds 
to the inverse correlation length of the statistical system. Usually masses are 
obtained from the exponential decay of suitable correlation functions. 

If a continuum limit with a finite physical mass exists, it means that by tuning 
the bare parameters (/5 in pure gauge theory) we can find a limit, where a-m 
goes to zero, while a physical quantity like a mass ratio ami/am2 remains finite. 
Thus the continuum limit of the lattice field theory corresponds to a critical point 
(second order phase transition) of a statistical system. When approaching the 
continuum, one has to increase L/a such that the correlators are not affected by 
finite-size effects. If one is interested in the continuum limit of explicit finite- 
volume quantities, L-m has to be kept fix as a goes to zero. 

In practice one picks one observable with small variance to set the scale, i.e. 
other observables are expressed in units of this one. For SU{?>) gauge theory and 
also QCD the most commonly used quantity is the hadronic length scale Tq [32], 
defined through the force F{r) between two static color sources. More precisely, 
ro is the distance where 



The choice of the constant 1.65 is based on phenomenological quark potential 
models, where the distance ro defined in this way corresponds to approximately 
0.5 fm. Assigning physical units to lattice results in this way introduces large 
systematic errors due to uncertainty in the physical value of tq, but as long as 
results are expressed in terms of tq this provides a well-defined way to compare 
lattice results from different actions and/or lattice spacings. For the simple pla- 
quette gauge action the value of r^/a is known quite precisely for the relevant 
range of (3 values [33]. 

While this amounts to a non-perturbative determination of a(/3), the evolu- 
tion a-dgo/da can also be calculated in perturbation theory. The result is valid for 
go going to zero and predicts an exponential decrease of the lattice spacing as go 
vanishes {asymptotic scaling). Thus the continuum limit is obtained by sending 
go to zero, or equivalently /3 to infinity. As will be detailed in the next section, 
physical observables approach their continuum values with powers of the lattice 
spacing a such that in practice a small range of (3 values is generally sufficient for 
a reliable continuum extrapolation. 




(3.13) 



3.3 Symanzik improvement 
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The situation in full QCD is more complicated since additional dimensionful 
parameters are present in the form of the quark masses. These have to be tuned 
such that in the (3 —>■ oo limit a corresponding number of renormalized quan- 
tities is kept constant. If a lattice regularization with (at least remnant) chiral 
symmetry is employed, a massless continuum limit is obtained by simply setting 
the bare quark masses to zero at each value of /5. In other cases (e.g. Wilson 
fermions, see Section additive mass renormalization requires that one finds a 
line in bare parameter space, where a suitable defined quark mass (i.e. through 
the PCAC relation) vanishes. 

3.3 Symanzik improvement 

A very insightful way of looking at the continuum limit of a lattice field theory is 
due to Symanzik [6-8]. Instead of taking the lattice as something that approxi- 
mates a continuum theory, he turns the tables and declares the discretized theory 
at finite cutoff to be the main object of interest. One can then construct a new 
theory in the continuum that asymptotically describes the lattice theory. 

More precisely, the discretization effects are modeled [34] by adding higher- 
dimensional interaction terms (accompanied by powers of the lattice spacing) to 
the original continuum action 



Here Cq is the continuum QCD Lagrangian and the terms £i,£2,--- are lin- 
ear combinations of those local operators, which The coefficients a, a^, ...are 
additional couplings with negative mass dimension, which renders the theory de- 
scribed by ()3.14|) non-renormalizable by power-counting. However, we are not 
interested in the renormalization of these "new couplings". Instead ()3.14|1 is 
treated as an effective theory to finite order in the latter (but to all orders in g^) 
and is thus renormalizable order by order in a. It is in this sense that renormal- 
ized correlation functions are matched between the lattice at finite cutoff and the 
effective continuum theory. 

This is analogous to the phenomenological approach of describing yet unde- 
tected substructures or the effects of heavy particles through higher-dimensional 
interaction terms. The most prominent example is probably Fermi's current- 
current interaction that approximates two weak-interaction vertices connected 
by the propagator of a. W boson. In the low-energy limit E ^ fn^, the W- 
propagator can be neglected and Fermi's description becomes exact. 

In this sense the additional terms in Symanzik's low-energy effective theory 
()3.14|) represent "new physics" (namely the lattice artifacts) entering at a scale 
of a~^. Naturally, this description can be valid only for energies small compared 
to the cutoff, i.e. lattice spacings fine compared to hadronic length scales. 




(3.14) 
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In order to study how correlation functions approach their continuum limit, 
one also has to represent renormalized (composite) lattice fields in the effective 
theory. These effective fields 

0cfr = 00 + a0i + a^02 + • • • (3.15) 

are given as a sum of terms of increasing mass dimension that have the same 
lattice symmetries as the (composite) lattice field to be represented. 

With the help of the effective continuum theory the connected renormalized 
n-point lattice correlation function 

Gn{xi, ...,Xn) = (Z^)"(</)(a;i) . . . </)(x„))con (3.16) 

can be written in terms of expectation values with the QCD Lagrangian only 

Gn{xi,. . . ,Xn) = (0o(a;i) . . . 0o(a;„))con 

4>o{Xn con 

n 

+a^(0i(xo) . ..(f)i{xk) ■ ■ ■(j)o{xn))con + O(a^) . (3.17) 
fc=l 

The higher-dimensional terms from the action and the composite field now appear 
as operator insertions, thus explicitly showing the origin of the lattice artifacts. 
If the symmetries of the lattice action in question are such that the terms Ci 
and 01 are forbidden, correlation functions will approach the continuum with 
a rate proportional to a^. This is the case for the lattice chiral symmetry of 
Ginsparg-Wilson fermions, which will be discussed in Section I^TTl 

The Symanzik improvement programme aims at a removal of the lattice ar- 
tifacts of a given action order by order in the lattice spacing. This is done by 
modifying the lattice action and composite fields such that the terms Ck and (pk 
in the effective continuum theory vanish up to a fixed order k. If e.g. Ci is a 
linear combination of operators Oi, improvement (of the action) to order a can 
be achieved by adding a term ^CjOj to the action, where Oi is a lattice repre- 
sentation of Oi and Cj are suitably chosen chosen improvement coefficients. The 
same procedure is also applied to the fields 0. 

While in principle the improvement programme can be applied to arbitrary 
orders in the lattice spacing, in practice the increasing number of possible im- 
provement terms limits its application to a cancellation of the cutoff effects linear 
in the lattice spacing. For Wilson fermions 0(a) improvement of the action and 
the isovector currents and densities will be discussed at the end of Section 14.21 



Chapter 4 
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In the path integral formulation fermions are represented by two Grassmann- 
valued fields, ip and ip. On the lattice these become an assignment of ip and ip 
to each lattice site = an^, n^j, G N^. They carry color . . .) and Dirac 

(A, B, . . .) indices. A lattice action S combines the two lattice vectors into a 
scalar and is hence defined in terms of a matrix M = D + rriQ, where D is the 
lattice Dirac operator and rriQ the bare quark mass 

S= Yl ^i^rAM{x,yrA%m'^. (4.1) 

x,y;a,l3;A,B 

In the following several choices for D will be discussed and in particular we will 
see which properties of the continuum Dirac operator jfj,D^ can be preserved on 
the lattice. 

4.1 Chiral symmetry on the lattice 

As we have seen in Chapter |21 the global chiral symmetry is a central aspect 
of continuum QCD. Unfortunately it also seems to be the symmetry that is 
most difficult to obtain in a lattice formulation. In their 1981 paper [35] Nielson 
and Ninomiya gave a mathematical foundation to this finding. If we denote the 
Fourier space representation of the translation-invariant Dirac operator D by 
D{p), the famous Nielson-Ninomiya no-go theorem states [26] that there is no 
lattice Dirac operator that simultaneously has all of the following properties. 

o [locality] D{p) is analytic in the momenta Pfj_ with period 27r/a, 

o [small p limit] For small momenta D{p) behaves like i'j^p^ + O(ap^), 

o [no doublers] D{p) is invertible at all non-zero momenta (mod 27r/a), 

o [chirality] D anti-commutes with 75. 



21 



22 



Lattice QCD 



The first condition ensures tliat D is an essentially local operator, i.e. that 
in position space the coupling between sites decreases exponentially with their 
distance. Obviously the small momentum and invertibility conditions are required 
to obtain a continuum limit with the right number of particles and the correct 
dispersion relation. Finally, the last condition guarantees the invariance under 
the chiral transformations discussed at the end of Section 12.11 The violations of 
these properties are useful to classify the various fermion discretizations. 

Directly transcribing the continuum action ()2.8|) on a space-time lattice using 
the covariant lattice derivatives ()3.3p leads to the naive lattice fermion action 

•S"™ = ^ V^(x)7^V/,^(x) + moilj{x)ip{x). (4.2) 

X 

This action is ultra-local and due to its Dirac structure it also inherited all the 
(isovector) chiral properties of its continuum counterpart. By construction, the 
behavior for small momenta is also correct. However, in the chiral limit mg 
the inverse momentum-space propagator 

S~^{p) =mo + -^'yf, sin(p^a). (4.3) 

of this action has zeros on all corners of the Brillouin zone in addition to the 
one ai p = (0, 0, 0, 0). These spurious poles of the propagator cannot be ignored, 
since they survive in the continuum limit. They describe additional unwanted 
particles, the so-called "doublers" and hence the naive lattice fermions violate 
the third property mentioned above. To make things worse, the doublers come 
in pairs with opposite chirality, thus spoiling the U{1)a axial anomaly [36].^ 

To achieve a correct continuum limit, Wilson added a second derivative term 
to the naive action ()4.2j) . This gives a mass of order to all doublers such 
that they decouple in the continuum. The locality of the action is not affected 
and also the small momentum behavior is only modified at O(ap^). The Wilson 
fermion action with Wilson parameter r is given by 

Sf = a'^y^^(x)(£)„ + mo)ip{x) 

X 

= ^ ip{x)-f^Vf,ij{x) - arij{x)\V ^*^^l){x) + mQi){x)ip{x) . (4.4) 

X 

The Wilson term proportional to r (0 < r < 1, [38]) is diagonal in Dirac space 
and hence violates the last point in the no-go theorem. The properties of Wilson 
fermions and in particular the consequences of this explicit breaking of chiral 

^One can define a different axiai transformation to remove tliis cancellation [37], but this 
still produces the anomaly for 16 flavors. 
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symmetry (also with respect to the axial anomaly) will be discussed in more 
detail in the next section. 

In 1982 Ginsparg and Wilson [39] suggested a way to evade the no-go theorem 
and preserve the consequences of chiral symmetry. They proposed to replace the 
condition {D, 75} = with the milder condition 

L)75 + = aD-f^D (GW relation)^ (4.5) 

An action defined with a Dirac operator satisfying this relation is invariant under 
a continuous chiral symmetry, which becomes the continuum chiral symmetry as 
the lattice spacing a goes to zero. Indeed, it is easily verified that ()4.5j) implies 
the invariance of i/jDip under the transformation 

= 75(1 - laD)ij and 6^; = tp{l - \aD)-i^ (4.6) 

and also its isovector counterpart. Comparing this to the axial transformation in 
()2.H()|1 . it can be seen as a redefinition of 75 to 75 = 75(1 — 0/^/2). 

Since one now has an exact chiral symmetry at finite lattice spacing, there 
exist exactly conserved (albeit rather complicated) isovector axial currents [40]. 
Also the axial anomaly with GW fermions appears not only in the continuum 
limit. Inserting 75 into ()2.3H) one immediately obtains the anomaly in the form 

{d,A,{x)0,^,) = -a(axtTr [75D]) . (4.7) 

In [26] it is shown that the GW relation also implies that for any gauge back- 
ground —all [75-D] = 2A^f X index(D), where A'f is the number of quark fiavors 
and the index is the difference of the right- and left-handed zero modes of D. The 
topological interpretation of the anomaly is then obtained via the Atiyah-Singer 
index theorem [41], which (for smooth gauge configurations) states that 

index(D) = ^ W^/ F;,F;. ■ (4-8) 

It took 15 years to find solutions to the GW relation for the interacting case of 
QCD. These now include the (classically perfect) fixed point action [42], domain 
wall fermions [43] and the related overlap formalism [44]. Tempted by the beauti- 
ful properties implied by the GW relation, attempts have been made at dynamical 
QCD simulations using overlap fermions [45-47]. However, currently available al- 
gorithms and computer resources essentially forbid a physical application except 
in the quenched case. 

While already the global chiral symmetry of QCD is a very important aspect 
of the theory, its role in the electroweak theory is even more fundamental. There 

^The original and more general form is actually D75 + 75I? — 2aDR'^^D, where R is some 
local operator. 
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the flavor (or weak isospin) SU (2) is promoted to a local gauge symmetry and its 
presence is therefore necessary to provide a renormalizable theory. The existence 
of solutions to the GW-relation is thus crucial to define a chiral gauge theory 
non-perturbatively (see [48] for a review). Hence they help in establishing the 
renormalizability of the gauge sector of the standard model beyond perturbation 
theory. 

The "staggered" or Kogut-Susskind (KS) fermions [49,50] represent a very 
different approach by trying to directly address the problem of fermion doubling. 
Through a unitary transformation of the fermionic fields the action is brought to 
a spin-diagonal form, which is used to reduce the doubling problem to a four- 
fold "taste" degeneracy. While KS fermions keep a remnant chiral symmetry, 
they do not completely solve the doubling problem and also the construction 
of operators with correct quantum numbers is far from trivial. An additional 
theoretical problem [51] with the locality of the staggered operator might lie 
in the "square-root-trick" [52], which is used to further reduce the number of 
fermion flavors. With respect to the anomaly, staggered fermions have the same 
properties as the naive ones, i.e. the species doubling spoils an anomaly for any 
number of flavors different from four (or multiples thereof). 

The last lattice fermion formulation we mention here is again based on the 
Wilson Dirac operator. The (in practice) most immediate consequence of its lack 
of chiral symmetry is the fact that it is not protected against zero-modes even 
at non-vanishing quark mass. In the so-called twisted mass formulation [53] of 
lattice QCD a term proportional to «75T'^ (with the Pauli matrix acting in 
flavor space) is added to the Wilson Dirac operator for two quark flavors. As a 
result this operator has a manifestly positive determinant. Also here part of the 
continuum chiral symmetry is recovered at the cost of breaking the vector flavor 
symmetry such that the total number of conserved currents remains the same as 
for (untwisted) Wilson fermions. 

One particularity of this lattice fermion formulation is related to a spurionic 
symmetry of the terms appearing at 0(a) in the Symanzik expansion ()3.17|1 . 
Using this symmetry it is possible to set up a twisted mass formulation such that 
0(a) lattice artifacts cancel in most physical observables without the need to 
tune improvement coefficients [54]. 

4.2 Wilson fermions 

We will now proceed with a more detailed discussion of Wilson's fermion action 
fl4.4|) . Inserting the expressions for the covariant derivatives fj3.3|) and fixing the 
Wilson parameter r to unity, it can be written as 




X 



(4.9) 



4.2 Wilson fermions 



25 



where a sum over fi is still implied in all but the last expression. The action can 
be rearranged into a diagonal and off-diagonal part 



+tp{x){mo + - )i;{x) } . (4.10) 



a 



In practice the hopping parameter k = (2amo + 8)~^ is used instead of the bare 
mass rriQ. If one rescales the fields ip and ^p hy a factor -\/2K/a^/^, the Wilson 
action assumes the simple form 

^/ =E { - E [^{x)il-l,)U,ix)i;{x+ft) + ij{x+ft){l+^,)Ul{x)i;{x) 

+%l){x)ip{x)^. (4.11) 

The transformations of charge conjugation C, parity P and time reversal T listed 
in Table 14.11 are all discrete symmetries of the Wilson fermion action (this is of 
course also true for the plaquette gauge action). In addition it also inherited the 
75-Hermiticity (751^75 = due to the use of the symmetric derivative V^) from 
the continuum. 

As already mentioned, chiral symmetry is broken at order a by the Wilson 
term and recovered only in the continuum. As a consequence the point of vanish- 
ing bare quark mass is no longer endowed with an enhanced symmetry (leading to 
the conservation of the axial current) and protected from renormalization. How- 
ever, by tuning the bare quark mass tuq to a non-zero value one can still find a 



c 


U,{x) ^ u;{x) 


i){x) ^ Cij^ix) 
i;{x) ^ -4j^{x)C-^ 
C = 27072 


p 

X Xp = (xo, -x) 


p 

Uo{x) — > Uo{xp) 
Uk{x) ^ Ul{xp-k) 


p 

ij{x) — > 1o^{xp) 

— p — 

il){x) — > i^{xp)-fo 


T 

X ^xt = (-a;o,x) 


Uoix) ^ U^oixT-0) 

T 

Uk{x) — > Uk{xp) 


tpi^x) TiP{xt) 
tplx) iI){xt)T~^ 
T = i7o75 



Table 4.1: C ,P and T transformations on the lattice. 
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point, where the axial current is conserved (up to cutoff effects). At this critical 
mass mc{gQ) a quark mass derived from the PCAC relation (see Section IK!T|) van- 
ishes and thus the theory can be parameterized in terms of the bare subtracted 
quark mass = mo — rric. In perturbation theory the critical mass if found to 
be [55] 

me(^7o') = -0.2701 (7o' - 0.0430 ^?o^ + 0{g^,) . (4.12) 

In addition to this additive quark mass renormalization the axial current itself 
requires a finite renormalization, which is discussed in Section 14.31 In analogy 
to Section [2.31 we now study the lattice currents associated with the vector and 
axial flavor transformations in more detail. 



4.2.1 Vector currents 

The isospin symmetry, which the continuum action possesses in the case of de- 
generate quarks, is not broken by the Wilson term. Hence a global isospin vector 
variation of the Wilson action vanishes and a local one, parameterized by u;"(x), 
can be written divergence 

5^Sf = a'Yl , (4.13) 

X 

where 9* is the backward lattice derivative and V^{x) the split-point vector 
current 

V^{x) = '^{^P{x+fi^T%l+^^)Ul{x)i^ix)-i^^^^^ . (4.14) 

Although one now has an exactly conserved vector current for Wilson fermions, 
in practice the local current 

V^{x) = ^(x)74rXx) (4.15) 

is generally used for convenience since sources for correlation functions of this 
current are easier to construct. Not being a Noether current, it is not protected 
from renormalization. However, its normalization factor Zylg^) is calculated es- 
sentially as a by-product of the axial current renormalization discussed in Chap- 
ter El 

4.2.2 Axial currents 

Of more interest are the axial transformations. Let us first consider the case 
of naive fermions, where also chiral symmetry is preserved on the lattice. The 
associated current A'^{x) is conserved at zero quark mass due to 

S^Sn = a^^w(x)[-9*I^(x)+2moP'^(a;) 

X 

with I^(x) = ||^/;(a; + /i)|r'^7^75f/i(x)^/'(x) + ^(x)|r''7^75?7^(x)V^(x+/i)} 
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At non-vanishing bare quark mass the divergence of this axial current is propor- 
tional to the pseudo-scalar density 

P"(x) = 7/^(x)75|rXx) . (4.16) 

In this way naive lattice fermions have the chiral symmetry of the continuum if 
the local axial current 

A';^{x) = V5(x)7^75irXa;) (4.17) 

is replaced by the split-point axial current A'J^{x). However, this symmetry is 
broken by the Wilson term such that even at vanishing bare quark mass no 
conserved axial current can be defined. Denoting the axial variation of the Wilson 
term by X°'{x) we obtain [56] 

S^Sf = a^Yl '^"^^) \ ~ + 2moP"(x) + X%x)\ , with 

X 

^""i^) = -^$^{^(a^)K75t/M(^)^(^ + A)+^(^ + A)K75f/;(x)V^(x) 

+[x^(x-/i) -4P^{x)}. (4.18) 

Here the explicit breaking of chiral symmetry by the Wilson term can be seen 
from the fact that X°- can not be written as the divergence of a current. 

If Wilson fermions are an admissible regularization of QCD, it must be possi- 
ble to reproduce the correct axial anomaly in the continuum limit. Performing an 
isosinglet axial transformation of the Wilson action one obtains eq. ()4.18j) with 
the replacement (|r° — > 1). It is precisely the isosinglet axial variation X of the 
Wilson term, which (under certain assumptions [57]) reproduces the correct axial 
anomaly 

lim X = e F"- (4 19) 

4.2.3 Improvement 

For the Wilson fermion action the Ci term in Symanzik's effective theory can 
contain the five expressions {cr^^ = ^['y^,jj^) 

O3 = mqTr {F^uF^,} , O4 = m^{4j^ ^^j - tPi^D^i;} , (4.20) 
O5 = mlipip , 

since these are compatible with the symmetries {C,P,T and 75-Hermiticity) of 
the Wilson-Dirac operator. Note that none of these terms can appear in the 
effective action for a Dirac operator that satisfies the GW-relation ()4.5|) . 

Through formal manipulations of the Euclidean functional integral (which 
amount to an application of the classical field equations) two of these (e.g. O2 
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and O4) can be eliminated [34]. Wliile this is necessary from a practical point of 
view, it restricts the improvement to on-shell quantities, i.e. correlation functions 
over finite distances. However, since in practice all physically relevant matrix 
elements are on-shell, this does not represent any loss of applicability 

On the other hand, adding and merely amounts to a rescaling of the 
bare coupling and quark mass by factors l+0(amq). We will return to this point 
in the discussion of normalization factors. What remains is the Pauli term Oi. 
On the lattice the field strength tensor can be represented by 

= ^[Q,u-Qu,], (4.21) 
where Qf,u{x) = U^-^^^ + + U^--^-^ + Ux■v-^, , (4.22) 



to arrive at the action [58] 



5'sw = ^'^(x) D^ + mo + ac^^la^^F^^ ip{x) . (4.23) 



The normalization of the improvement term is chosen such that the improvement 
coefficient Csw = 1 at tree-level in perturbation theory. In general however, 
Csw depends on and a non-perturbative determination for the two-flavor case 
[59, 60] is summarized by the interpolating formula 

_ 1 - 0.454^g - 0.175^0^ + 0.012^6 ^ 045^8 
CsAQo)- 1 - 0.720(72 • ^^■^^> 

For the gauge action Ci vanishes and thus at 0(a) no improvement is needed. 

The number of possible improvement terms for the isovector axial and vector 
currents can be reduced using the same arguments as for the Wilson action. The 
on-shell improved currents are given by 

{AT, = A^ + acA9^P% (4.25) 

(^i)^ = V^ + acydXu, (4.26) 
where T;^(x) = #(x)a,,4rXx) . (4.27) 

The pseudo-scalar density P"'{x) is already improved, i.e. no 0(a) counter term 
with the correct symmetries exists. For the vector current correlation functions 
we will consider, the term proportional to cy does not contribute. The 1-loop 
result for ca is [61] 

ca(^7o) = -0.00756(7o' + 0{g^) (4.28) 



and a non-perturbative determination will be presented in Chapter 
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4.3 Current renormalization 

A scheme, in which the renormahzation conditions are imposed on a set of correla- 
tion functions at a point {go, arrio) in bare parameter space, would result in renor- 
malization factors, which depend on both the bare coupling and the bare mass 
as well as a renormalization scale a/x. Inherently simpler are mass-independent 
schemes, where one imposes the normalization condition at zero quark mass with 
the consequence that no dependence of the renormalization factors on the quark 
mass is introduced [62]. 

Neglecting 0(a) improvement, in such a scheme one would introduce renor- 
malized parameters through e.g. mR = Zmfriq, where the subtracted mass 
mq = rriQ — rric parameterizes deviations from the critical line and the renor- 
malization factor depends on the bare coupling Qq as well as the renormalization 
scale a/i. As is discussed in [63], such a definition leads to uncanceled 0(am) 
effects and is thus not compatible with the improvement programme. This can 
be understood from the fact that in Section 14.2.31 we ignored 0(a) counterterms, 
which merely amounted to a rescaling of the bare parameters. 

If we want to employ a massless renormalization scheme and at the same 
time have on-shell 0(a) improvement, we need to allow for a more general form 
of renormalized quantities. These are expressed in terms of the modified bare 
coupling and quark mass 

= 9o{^ + ham^) , (4.29) 
rriq = "^q(l + ^momq) , (4.30) 

where the 6-coefficients are functions of Qq. The renormalized coupling and quark 
mass can then be written as 

= 9lZg{gl,a^x) (4.31) 
mR = m^Zm{gl,a^i) . (4.32) 

Also for local composite fields the l + 0(amq) counterterm is conventionally not 
included in the definition of the improved field 0i, but rather in the renormalized 
field according to 

0r(x) = Z^{gl, a/i)(l + 6^amq)0i(x) . (4.33) 

In this context the isovector axial and vector currents are special since in the 
(massless) continuum limit they become the Noether currents of the flavor chiral 
symmetry and at finite lattice spacing their normalization can thus be fixed by 
imposing a continuum Ward identity. The latter is a local identity and hence the 
renormalization factors of the isovector axial and vector currents do not depend 
on a scale. In addition they approach unity for going to zero, i.e. when the 
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flavor chiral symmetry is recovered [56,64]. 



{A^r, = ZA{fo)a + bAam^){Ar,, (4.34) 

{VnT, = Zy{~gl){l + byam^){V,r^. (4.35) 

To one loop in perturbation theory [65-67] one finds tliat 

Za = 1 - 0.116458 (?o' + ((7^) , (4.36) 

Zy = 1 - 0.129430 gl + 0(^^) . (4.37) 



One should point out again that we have ignored the existence of an exactly 
conserved vector current ()4.14|1 for Wilson fermions and treat and on 
exactly the footing. The normalization of the split-point vector current ()4. 14j) is 
given naturally, i.e. it does not renormalize. 

The renormalization of the pseudo-scalar density 

{PnTix) = Zp(l + 6pam,)P'^(x) (4.38) 

on the other hand will depend on the scale, at which the normalization condition 
is imposed, i.e. Zp = Zp(^Q,a/i). 

A massless scheme requires that the normalization conditions are defined at 
vanishing quark mass. In a non-perturbative implementation this implies that 
numerical simulations have to be performed at zero (or at least very small) quark 
mass. With periodic boundary conditions such simulations are technically very 
difficult or even impossible since there is no lower bound on the spectrum of 
the Dirac operator. This problem is cured in the Schrodinger functional setup, 
which will be discussed in the following chapter, at the end of which we will 
define a renormalized quark mass in terms of correlation functions of renormalized 
composite fields. 



Chapter 5 

The Schrodinger functional 



So far no boundary conditions have been specified since tliey would have been of 
no importance in the discussion of the local symmetries. It was thus implicitly 
assumed that we are working either in infinite volume or in finite volume with 
periodic boundary conditions for the fields. All numerical results report here 
were obtained in the Schrodinger functional (SF) setup, which employs periodic 
boundary conditions in the three spatial directions and fixed (Dirichlet) boundary 
conditions in time. 

There are several advantages of this method compared to the torus, where 
all four directions are periodic. The Dirichlet boundary conditions provide an 
infrared cutoff (inversely proportional to the time extension T) to the Dirac 
operator [10].^ The gluonic boundary conditions can be used to induce a color 
background field, which was employed e.g. in the definition of a running coupling 
[2,68]. The fermionic boundary fields can serve as sources for correlation functions 
in an elegant way, particularly if a zero-momentum projection or wave-functions 
are required. This results in mesonic correlation functions with smaller statistical 
error than on the torus. 

A vast literature (see e.g. [9,10,34,69]) exists on the SF and we will therefore 
restrict the discussion here to the aspects required in the following chapters, in 
particular the ingredients required for the definition of the lattice SF with Wilson 
fermions. Some SF correlation functions needed in the axial current improvement 
(Chapter [7j) and renormalization (Chapter |H)) are also introduced. 

The QCD Schrodinger functional is the propagation kernel of a field configu- 
ration C, p, p at Euclidean time to a configuration C\p\ p' at time T. Following 
Feynman, this is written as an integration over all interpolating field configura- 
tions 





This can be rigorously shown only in the ^ Q limit. 
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5.1 Lattice SF with Wilson fermions 

We consider an L'^ x T space-time cylinder, i.e. the set of sites {k = 1, 2, 3) 

{x I x/a e Z^0 < xo < r,0 < Xfc < L} (5.2) 

and the associated links [/^(x).^ The gauge links are spatially strictly periodic, 
while the fermionic fields are allowed to pick up a phase 6k at the boundary 

U^{x + Lk/a) = U^{x) 

-^{x + Lk/a) = e'^'ipix) (5.3) 
^p{x + Lk/a) = e-'^'^-ipix) . 

The Dirichlet boundary conditions are implemented via 

f^fc(a;) 1x0=0 = exp(Cfc) and [/fc(x)|^o=r = exp(C^) , (5.4) 

where Ck is a traceless anti-Hermitian 3x3 matrix. Throughout this work we 
use Ck = CI = 0, i.e. no background field. From the lattice action it will 
be obvious that only half of the fermion fields components on the boundaries 
actually couple to the interior of the SF cylinder. The SF therefore becomes a 
well-defined boundary value problem by prescribing {P± = (1 ± 7o)/2) 

P+^(x)|^^^o = p(x) , 7/;(x)P_|^.^^„ = p(x) , 
P_V^(x)|^^^^ = p'(x), ^(x)P+|^^^^ = p'(x). 

While the fermionic boundary fields (p, . . .) are always set to zero, functional 
derivatives of the effective action with respect to them are used to define correla- 
tion functions involving the SF boundary. Thus, expectation values of composite 
fields O in the lattice Schrodinger functional are obtained from 

{O) = \ ^ [vUV^Vi: O e-'^^^'^'^^] . (5.6) 



Z 



p'=p'=p=p=0 



In addition to the dynamical variables U, ip and ip, the composite fields may then 
also involve the "boundary fields" (see Appendix |X] for details) 

C(x) = , C(x) ^ 



^P(x) 5p(x) 

When acting on the Boltzmann factor in ()5.6j) . these have the effect of insert- 
ing certain combinations of and tp close to the boundary, together with the 

^Except for Uo{x) when Xq = T. 
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appropriate gauge field variables to ensure gauge covariance. While the lattice 
action S[U, ijj, ijj] is in principle given by ()H.4|) and ()4.2H|1 together with the above 
boundary conditions, the issue of 0(a) improvement needs to be re-addressed 
due to boundary effects. 

In the presence of the Dirichlet boundary conditions, surface terms have to 
be added to the effective action ()3.14p . At 0(a) this is [63] 

lim J d3x{i3i(x)U„=, + B[{^)\,,=T-e} , (5.8) 

where Bi and B[ are linear combinations of composite fields of mass dimension 
four. For the gauge action the only fields that can contribute are Tr {F^iFki} and 
Tr{Fofc-Fofc}. On the lattice these are implemented by modifying the weight of 
the plaquettes in the boundaries (cs/2) and those touching the boundaries (ct). 

On-shell 0(a) improvement of the fermion action can be achieved by e.g. 
adding the field products 



at Xq = and 



Oi2 = i'P-D.i, + i^DoP+iJ , 

Ois = ^P+ikDkiJ - ijDkikP-ij , 



at xq = T. The corresponding lattice discretizations of these terms are added to 
the action with coefficients Ct for the temporal derivatives and Cs for the spatial 
derivatives. 

To keep the discussion transparent, the explicit form of the improved action 
is given in Appendix ^ Here we want to note only that for our specific choice of 
boundary conditions the terms proportional to Cg and Cg do not contribute, while 
Ct and Ct are known perturbatively [70,71] 

Ct = 1 -0.089^^-0.03^^, (5.11) 
Ct = 1- 0.018 (?o'. (5.12) 

Unless T/a is very small, a non-perturbative knowledge of these coefficients is 
not necessary since the lattice artifacts they modify are surface effects and thus 
suppressed by a/T compared to the cutoff-effects from the bulk. Moreover, in [63] 
it is shown explicitly that these term can not contribute 0(a) terms to matrix 
elements of the PCAC relation (since the latter holds exactly in the continuum 
theory). 

The two-point functions are derived by taking functional derivatives of the 
generating functional and all the fermionic propagators involving the boundary 
and bulk quark fields are also given in Appendix 1X1 
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5.2 Fermionic correlation functions 

We are interested in mesonic correlation functions, i.e. we want to excite states 

with the quantum numbers of isovector pseudo-scalar particles. In terms of the 
boundary fields ()5.7|) possible choices are 

O'' = a6^C(u)75|r'^C(v), (5.13) 

u,v 

(9- = a6 5^C'(u)75|rr(v). (5.14) 



u,v 



With these operators both constituent quarks are projected to zero momentum 
separately. In Chapter [7| we will also consider generalizations of ()5.131 I5.14|) . 
which make use of spatial trial wave functions to infiuence the excitations of 
pseudo-scalar states. 

To obtain non-trivial expectation values in ()5.(ij) the operator appearing there 
needs to be compatible with the lattice symmetries, in particular it must be invari- 
ant under parity and an isoscalar. Restricting ourselves to two quark bilinears, 
there are five possible combinations involving ()5.13l I5.14|) 

/A(a;o) -- 

fpixo) -- 

gA{xo) -- 

9p{xo) -- 

and /i -- 



X 


(5.15) 


X 


(5.16) 


X 


(5.17) 


X 


(5.18) 




(5.19) 



Their explicit form in terms of quark propagators is given in Appendix lB.il From 
these correlation functions one can define a current quark mass m through the 
PCAC relation ()2.28p . If the transformations ()2.17|) are applied in the interior, 
the boundary fields O and O' are "external" and we get 

Bf^iA'^^O'') = 2m(P°C") (5.20) 

and thus e.g. 

2/p(T/2) ^'-'^^ 

defines a bare current quark mass since the spatial derivatives vanish under pe- 
riodic boundary conditions. The lattice artifacts, by which eq. ()5.21|1 is affected 
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can be reduced by using the improved axial current ()4.25|1 to arrive at 

go/A(T/2) + acAaogoVp(T/2) 

2/p(T/2) ^^-''^ 

as a possible estimator for an improved current quark mass. Of course, also qa and 
gp could be used instead. In fact, for vanishing background field {Ck = C'f, = 0) 
time-refiection is a symmetry of the SF and fx{xo) = gx(T—xo). An average of 
/ and g can then be used to improve the statistical precision. 

Since their derivation is based on local symmetry transformations of the con- 
tinuum action, the quark masses defined in this way are independent of the lattice 
extensions L and T up to cutoff effects of O(a^) in the improved theory. This is 
explicitly shown in [61]. 

The Schrodinger functional also defines a (finite volume) renormalization 
scheme, which has been used e.g. in the calculation of the scale dependence 
of the SF coupling g'^ [2]. In the SF scheme the renormalization scale fi is iden- 
tified with the inverse box length L~^. In the context of a project to determine 
the scale evolution of quark masses, the renormalization constant of the pseudo- 
scalar density has been determined non-perturbatively in this scheme [13, 14]. 
To this end one requires that a certain matrix element of (-Pr)" is equal to its 
tree-level value c. The concrete requirement is that 

Zp{go,L/a) = ^j^^j^ (5.23) 

at vanishing quark mass. One can then proceed and insert the renormalized 
currents into ()5.22|) to define a renormalized improved quark mass 

_ ZA(l + feAamq) 2n /c.o^^ 

■"^ = -m-. — ; mi + 0(a) (5.24) 

Zp(l + 6pamq) ^ ' ^ ' 

= ^{l + {bA-bp)am^)mj + 0{a^) . (5.25) 

Zp 

Since depends only on the (modified) bare coupling, the entire scale (and 
scheme) dependence of the renormalized quark mass enters through the normal- 
ization of the pseudo-scalar density. 

The evolution of m is then traced to very high energies, where the SF scheme 
can be related to a more common renormalization scheme like MS by means of 
a 1-loop calculation. The application of this programme to the strange quark 
mass [15] is one of the immediate applications of the axial current improvement 
and renormalization presented here. In the quenched approximation this has 
been done in [13,72] 



Chapter 6 
Algorithmic issues 



The basics of our lattice theory as well as the framework, in which we implement 
the axial current improvement and renormalization, have now been laid out. 
However, before proceeding to a detailed discussion of the latter, it is necessary 
to detour into the more technical aspects of full QCD simulations. 



6.1 Hybrid Monte Carlo algorithms 

Since the fermionic action is bilinear in the quark fields, the Grassmann integral 
can be performed analytically, resulting in the determinant of the Dirac operator. 
This step is even necessary unless one can find an efficient way of dealing with the 
Grassmann fields directly. For two quark fiavors one has (writing M = D + m) 

Z = lvUVi)V^ exp [-Sg[U]-^ i)M[U]^'^ (6.1) 

flavors 

= jvU e~^«t^l det M[Uf . (6.2) 

Due to the 75-Hermiticity of M this can be expressed in terms of the positive 
matrix = Q^Q with Q = 75M, such that detM^ = det QQ^. With the intro- 
duction of a so-called pseudo-fermion field the determinant can be expressed 
bosonic Gaussian integral 

Z = jvUVcj) exp {-Sg[U] - 0^(QQ^)"V} = jvUVcj) e'^^^ . (6.3) 

To generate gauge field configurations distributed according to the effective ac- 
tion, the hybrid Monte Carlo (HMC) [73] algorithm evolves the gauge field in 
a fictitious time r (molecular dynamics). This evolution is described by anti- 
Hermitian traceless momentum matrices vr^(a;) 



d_ 



U^{x) = ■n^{x)U^{x) . (6.4) 
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Adding a Gaussian action for tlie momenta to Scs, we obtain tlie Hamiltonian of 
the molecular dynamics (MD) evolution 



The evolution equation of the momenta vr is derived by requiring the Hamiltonian 
to be a constant of the MD motion. The update of the momenta requires the 
computation of the "force", i.e. 6Sca/SU, which involves an inversion of Q^Q for 
the fermionic contribution. For this inversion the conjugate gradient algorithm is 
used (see Section lU^ . Since the momenta don't couple to the gauge or pseudo- 
fermion fields they will have no effects on correlation functions involving U and 
(j). An HMC update step (trajectory) now proceeds from a starting configuration 
U{0) according to 

1. Generate a complex Gaussian random vector and apply Q to it, 
to obtain (p distributed according to 

2. Generate momenta according to ()6.5|) by drawing Gaussian com- 
ponents in the Lie algebra generators; 

3. For fixed 0, numerically evolve U and vr for a MD time tq using 
a numerical approximation of Hamilton's equations of motion; 

4. Accept U{to) with a probability e"^^ , where 



In case of rejection U{0) becomes also the next configuration in the Markov 
chain. In the next trajectory new pseudo-fermion fields are used to estimate the 
determinant contribution and also the direction of evolution is again chosen ran- 
domly by using new Gaussian momenta {refreshed molecular dynamics). Several 
aspects of the HMC are worth mentioning at this points. 

For the algorithm to satisfy detailed balance with respect to the effective ac- 
tion, the numerical integration scheme needs to be reversible and area preserving 
(in phase space). The molecular dynamics evolution is the computationally most 
expensive part of the HMC algorithm. 

The simplest reversible and area preserving integration scheme, is the leapfrog. 
With 6t = To/ustep this is defined through 



i7 = S'eff-iTr7r2. 



(6.5) 



AH = H{to) - H{0) (Metropolis step). 



71 



7v-6S^s/SU-6t/2 




(6.6) 



U 



TT 



exp{7i6r)U 

7T -6S^s/6U-6r/2 . 



Throughout this work we always consider molecular dynamics trajectories of unit 
length, i.e. Tq = 1 and 6t = l/ngtop- The half-steps differ from the exact integra- 
tion of Hamilton's equation by errors of order (5r^, whereas the nstep = 0(l/(5r) 
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Figure 6.1: Left: Average plaquette from different algorithms in SU{2) gauge 
theory on a 6^ lattice with periodic boundary conditions at /3 = 3.5. Right: scaling 
of the error with CPU time icpu for heatbath and HMC in this setup. 



intermediate steps each have errors of 6t^. More elaborate (and hence more ex- 
pensive) schemes can increase the acceptance rate by lowering the violation of 
energy (Hamiltonian) conservation. If the accept/reject step at the end of the 
HMC trajectory is omitted, one generates a systematically wrong ensemble and 
an extrapolation to 6t = becomes necessary (0-algorithm [74]). In this case 
observables are expected to differ by 0(5r^) from their correct value if a leapfrog 
integrator is employed. 

While most of the computer time is spent in evaluating the fermionic force, 
one should note that already for pure gauge theory, the HMC itself is a rather 
inefficient algorithm. Although it globally moves the field configurations, in prac- 
tice the autocorrelation times are large. 

To illustrate these points, results from three-dimensional SU (2) gauge theory 
are shown in Fig. 16.11 At this statistical precision the 0(5r^) errors from the 
leapfrog integration are clearly visible for the 0-algorithm (left plot). The HMC 
results for the average plaquette agree with the (much more precise) heatbath 
data for all values of the leapfrog step-size. A comparison of the statistical 
error as a function of the CPU time (right plot) shows that the heatbath gives a 
significantly more precise result for a given computational effort. 

Many improvements to the "basic" HMC described here have been proposed 
and tested to speed up full QCD simulations. This includes preconditioning and 
different pseudo-fermion representations of the fermion matrix as well as higher- 
order integration schemes for the MD evolution. However, so far none of these 
represent a real breakthrough in algorithmic development since the performance 
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gain is in general quite moderate (i.e. factors of < 2 at most). The situation 
might be different for a recent promising proposal [75], which relies on a domain- 
decomposition of the fermion matrix. Here significant advantages in the region 
of small quark masses are expected. 

Another possibility is to use in the MD evolution a different action (guidance 
Hamiltonian) than one really wants to simulate. This can be corrected either 
through an adjustment of the acceptance probability in the Metropolis step or by 
calculating a reweighting factor to be included in the data analysis. The latter 
possibility is employed in the polynomial hybrid Monte Carlo (PHMC) algorithm, 
which will be discussed in Section 16.61 

All of the simulations presented here are done using even-odd preconditioning 
[76], where the lattice sites are divided into an even and an odd part according 
to Xl/^^M- Since the hopping term in the Wilson Dirac operator couples only 
even to odd sites, it is possible to rewrite Q in a form, where the non-trivial 
part acts only on the odd sites. As a result the effective condition number is 
reduced and inversions become computationally cheaper. We will denote the 
even-odd-preconditioned Hermitian Dirac operator by Q. 

Another improvement is the introduction of a second pseudo-fermion field 
as proposed in [77]. More precisely, the implementation we employed is the one 
discussed in [78], where the pseudo-fermion action from ()6.3|) is split according 
to 

Sf. = 0l(Q^ + p^)-Vi (6.7) 
and Sf, = 4iP~^ + Q~^)4>2 , (6.8) 

such that 

detg2 = yP(/)iP(/)2 e-^^i-^^2 . (6.9) 

With this split-up the fermionic forces in the MD evolution become smaller, such 
that larger step-sizes can be used in the numerical integration. 

Finally, an improved Sexton- Weingarten [79] integrator is used in the molecu- 
lar dynamics. It uses different time scales for the gauge and the fermionic part of 
the forces. A performance improvement results since the computationally cheaper 
gauge forces are evaluated more often. This scheme partially removes the inte- 
gration errors of 0{6t^). In practice it is found that the Metropolis acceptance 
rate Pace behaves as 

-logPacc^l-PaccCX^rS (6.10) 

thus indicating that the errors of 0(5r'^) dominate. 

An important quantity to monitor in an HMC simulation is the Hamiltonian 
violation AH during one molecular dynamics trajectory. In the large volume 
limit its distribution is Gaussian with mean and width related through [80] 



2{AH) = (AH^) - {AHf 



(6.11) 
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This relation holds for the simulation described in Table m| while for the simu- 
lations in smaller volume (Tables IfTD and iD.ljl deviations were observed. Only in 
the limit of an exact integration, 6t ^ 0, the mean (and hence also the variance) 
of AH vanish. However, from the detailed balance condition it is possible to 
show that even at finite step-size [81] 

(e-^^) = 1 (6.12) 

on the generated ensemble (including of course the rejected proposals). A viola- 
tion of this condition could indicate problems with reversibility in the numerical 
integration of the equations of motion. 



6.2 Inverting the Dirac operator 

To evaluate the fermionic contribution to the MD force and to calculate fermionic 
correlation functions, matrix elements of the quark propagator S = need to 
be calculated on the gauge background produced in the Monte Carlo run. De- 
pending on the specific correlator we need different linear combinations of matrix 
elements. These are obtained by inverting M on different sources q. Writing the 
Dirac indices A,B,... and color indices a, (3, . . . explicitly, we need to solve 

M{x, y)l% ■ S{y, zf^l = g(x, z)Z • (6-13) 

Since the source for the inversion can only be a vector, we need one inversion 
for each combination of the "column" indices z, 7 and C of q. Exploiting the 
linearity of ()6.13p . one can add different sources to obtain the corresponding 
linear combination of columns of S with a single inversion. For a given source y 
the solution x of Mx = y is found by minimizing the residue 



\\Mx - 


-y\ 


12 




\y\ 


2 



(6.14) 



Since M is a large sparse matrix, iterative conjugate gradient methods provide 
an efficient tool for this task. The most commonly used version is the stabilized 
biconjugate gradient (BiCGstab) [82]. However, we found that in cases where 
the matrix M is rather ill-conditioned, it can happen that the BiCGstab does 
not converge to a solution. 

This poses a problem in particular because during the iteration, the residue 
r is not computed according to ()6.14j) . To save an additional application of 
M, the current residue is computed from its value during the previous iteration 
(iterated residue). Thus, due to accumulation of roundoff the iterated residue 
might become small, indicating that a solution has been found, while the real 
residue is still of 0(1). 
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Figure 6.2: Solver residues as a function of the iteration number. For difficult 
inversions (lower plot) the BiCGstab might indicate a convergence in the iterated 
residue although no solution has been found. 



The simple conjugate gradient (CGNE)^ solver has a much smoother conver- 
gence behavior and is very stable even when applied to ill-conditioned problems. 
Fig. 16.21 shows the real and the iterated residue for BiCGstab and CGNE during 
the iterations of a normal inversion and one, where BiCGstab did not converge 
to a solution. For CGNE the real and iterated residues start to differ only when 
the single-precision limit is approached.^ 

To account for the known insufficiencies of the BiCGstab while still making 
use of its superior performance, we perform the inversions for the calculation of 
correlation functions in the following way. We start by first running the BiCGstab 
solver until the required precision or the maximum number of iterations, O(IOO), 
is reached. This is followed by a CGNE inversion starting from the BiCGstab 
result with a maximum number of 0(1000) iterations, although for a "normal" 
inversion only 0(1) CGNE iterations are performed until the required precision 
is obtained. To compensate for precision loss due to roundoff, the CGNE is 
restarted up to ten times with the previously found solution. The real residue is 
always recomputed after the entire inversion and monitored in the data analysis 
program. 



^NE stands for normal equation, AI'^Mx = M^y. CG itself applies only to positive matrices. 
^Apart from global sums all our calculations are carried out in single-precision arithmetics. 
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6.3 Data analysis 

As mentioned earlier, the field configurations in a Markov chain produced by a 
hybrid Monte Carlo algorithm can be strongly (auto) correlated in Monte Carlo 
time. This needs to be taken into account in order to obtain reliable error esti- 
mates for quantities measured in such simulations. 

Traditionally this is done by combining a number of consecutive measurements 
into bins and then performing a naive error analysis on the binned data (jackknife 
analysis). However, it has been shown [83] that this procedure is only suboptimal 
concerning the reliability of the error estimate. 

For all our data analysis we use the explicit numerical integration of the 
autocorrelation function for primary and derived observables from [83]. The 
integrated autocorrelation time Tint enters the error estimate a in the form that 
the number of estimates N is "effectively" reduced according to 

= , (6.15) 

where u is the variance of the quantity in question. While the latter is a prop- 
erty of the observable and hence of the theory one is simulating, the size of the 
statistical error is also influenced by Tint, which depends on the algorithm used 
to generate the Markov chain. This method also provides us with an error of 
Tint and hence we can directly assess the reliability of our error estimate, i.e. the 
error of the error. 



6.4 Sampling problems on coarse lattices 

In order to motivate the choice of simulation algorithms in the evaluation of the 
axial current normalization condition (Chapter |HJ we now proceed to a more 
detailed discussion of the algorithmic difficulties we faced. This concerns the 
simulations at the coarsest lattice spacing we consider, a~0.1 fm at a bare gauge 
coupling of P = 5.2. 

It has long been established that in the corresponding quenched situation 
cutoff-effects are compatible with Symanzik's description (Section 13. 3j) and a 
continuum extrapolation can be started there. However, over the last years more 
and more evidence has been accumulated that for dynamical improved Wilson 
fermions at this lattice spacing the cutoff-effects are much larger than expected. 
As an extreme example, for three flavors the existence of a phase transition in 
the /J-K-plane has been numerically conjectured and is interpreted as a lattice 
artifact [84]. A summary of large scaling violations in the two-flavor-theory is 
given in ref. [85]. 

In the following we will establish that the algorithmic problems we encoun- 
tered in the MD integration and the efficient simulation of the canonical ensemble 
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can also be interpreted as cutoff-effects. More precisely, at the infrared end of 
the spectrum of the Dirac operator we find a behavior different from what one 
would expect close to the continuum. 

One of the important results will be our finding that the distortion of the Dirac 
spectrum makes it advantageous to deviate from importance sampling. In this 
context we study the behavior of the polynomial hybrid Monte Carlo (PHMC) 
algorithm [86, 87] in this situation and find it a very useful tool for a detailed 
investigation of the properties of the small eigenvalues of the Dirac operator. 

In Table 16.11 we list the lattice sizes and bare parameters for the simulations 
discussed in this chapter. In all cases we have T = 9/4 L for the time extension T. 
The bare quark mass m is defined through ()5.22|1 . where we used the perturbative 
value of ca from eq. ()4.28j) . In the algorithm column 'H2' refers to the HMC with 
two pseudo-fermion fields described at the end of Section 16.11 and 'Pn' stands 
for PHMC with a polynomial of degree n. This algorithm will be introduced 
in Section 16.61 To gain more statistics, each simulation was performed with 16 
independent replica. 



run 


L/a 


/? 




Lm 


algo. 


A'traj St 


p 

^ acc 


[Al] 


8 


5.2 


0.13550 


0.205(10) 


H2 


16-500 1/16 


91% 


[A2] 


8 


5.2 


0.13515 


0.307(9) 


H2 


16-520 1/25 


97% 


[A3] 


8 


5.2 


0.13515 


0.314(8) 


Pl40 


16-500 1/26 


87% 


[A4] 


8 


5.2 


0.13550 


0.195(7) 


Pl40 


16-400 1/25 


85% 


[A5] 


8 


6.0 


0.13421 


0.193(3) 




— quenched — 




[A6] 


12 


5.5 


0.13606 


0.287(3) 


H2 


16-240 1/20 


91% 


[A7] 


12 


6.26 


0.13495 


0.295(3) 




— quenched — 





Table 6.1: Summary of simulation parameters. 



6.5 Instabilities in the molecular dynamics 
integration 

We have seen that in the numerical approximation to the molecular dynamics evo- 
lution the Hamiltonian is conserved up to powers of the step-size 6t employed 
in the integration. Apart from these (usually small) deviations, under certain 
conditions the currently used integration schemes can become unstable and pro- 
duce very large Hamiltonian violations AH. For a more detailed discussion see 
ref . [88] , where a connection between these instabilities and large driving forces in 
the MD is proposed in analogy to a harmonic oscillator model. In this model the 
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5 10 45 50 

max. ^rejections in a row 



Figure 6.3: from run [Al]. Left plot: histogram. Right plot: our estimates 

of Tint in units of MD time separately for the 16 replica. In our normalization 
■^int = 0.5 means no autocorrelation (dotted line). 



integrator becomes unstable when the product of the force and the integration 
step-size exceeds a certain threshold. 

Large (positive) values of AH result in a rejection of the configuration pro- 
posed by the MD. In a histogram of e~^^ these contribute to bins close to zero 
while the distribution is peaked around one, indicating that in most cases the 
numerical integrator performs as expected. They can also lead to an unusual 
autocorrelation of this quantity, making the Monte Carlo error estimate difficult. 
In particular this applies also to the integrated autocorrelation time of e~^^ it- 
self. This is due to the long periods of rejection in the Metropolis step, which 
sometimes follow large AH values. 

Fig. Ifi.HI shows a histogram of e~^^ and also its integrated autocorrelation 
time from one of our simulations. In this data set there are several series of 
large AH values, during which the proposed configurations were rejected. In the 
distribution of e~^^ these lead to an additional peak close to zero. One also sees 
from the right-hand plot that e~^^ is noticeably autocorrelated only when a 
large number of proposals were rejected in a row. As argued above in these cases 
the error of Tint could be underestimated. These two effects might cause some 
concerns when using (e~^^) — 1 as an indicator for the absence of reversibility 
violations [78]. 

Spikes in AH have been observed by several collaborations using (improved) 
Wilson fermions in various setups (e.g. different gauge actions and volumes) at 
relatively large lattice spacings [59,88-90]. There these spikes have been traced 
back to large values of the driving force in the MD evolution and also their 
dependence on the quark mass has been investigated. 

Here we want to clarify a point, which is essentially implied by the previous 
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Figure 6.4: Monte Carlo history for one replicum of run [Al] with a long period of 
rejection. Configurations where the new proposal was rejected are marked by a dot. 
At r~ 340 the algorithm gets stuck with a configuration carrying an exceptionally 
low smallest eigenvalue Amin of Q^- 



observations [89,91], namely the strong correlation between spikes in Ai7 and 
small eigenvalues of the Dirac operator.^ In this way we hope to be able to 
separate physical effects from cutoff-effects, i.e. the occurrence of unphysically 
small eigenvalues. The low-lying eigenvalues are computed using the method de- 
scribed in [92]. In Fig. 16. 41 we clearly see a long period of rejection (corresponding 
to the rightmost data point in Fig. caused by the presence of a very small 
eigenvalue. Although we did not measure them, this is expected to produce large 
fermionic contributions to the driving forces since they involve an inverse power 
of the Dirac operator. 

We found the observed average Amin to be close to its tree-level estimate with 
Schrodinger functional boundary conditions [10]. However, the smallest Amin is 
an order of magnitude below that and we therefore consider these eigenvalues 
unphysical and will later establish their nature as cutoff-effects. 

Finally, following the procedure of ref. [78], the absence of global reversibil- 
ity violations is explicitly verified even for trajectories resulting in large values of 
AH. Nevertheless our experience shows that the increased cost of using a smaller 
6t such that no long periods of rejection occur is more than compensated by the 
reduction in autocorrelation time of all observables. The reason is that already a 
small decrease of the integration step-size greatly reduces the Hamiltonian viola- 

■^Here and in the following we will always refer to the eigenvalues of the square of the 
Hermitian even-odd preconditioned Dirac operator in the Schrodinger functional. For its 
precise definition see ref. [78] . 
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T 

Figure 6.5: Monte Carlo history for the A'rep = 16 rephca of run [A2] showing the 
normahzed [i;^i]f- 



tions. For example, repeating run [Al] with a step-size of 1/20 instead of 1/16, 
the longest period of rejection was 4 (instead of 47) consecutive trajectories. Such 
a rapidly changing behavior also supports the picture of integrator instabilities 
from [88] , where an abrupt increase of Hamiltonian violations is predicted when 
a certain 6t threshold is exceeded. 

6.6 MC estimates of fermionic observables 

We concluded in the previous section that unphysically small eigenvalues of 
produce algorithmic problems only on a practical and not on a theoretical level. 
But apart from slowing down the algorithm these small eigenvalues also cause 
problems in the MC evaluation of fermionic Green's functions. 

As an example we consider the Schrodinger functional correlation function /i 
as defined in eq. ()5.19|) . It is the correlation between pseudo-scalar composite 
fields at the first and last time-slice, respectively. We will denote its value on a 
given gauge field configuration by [0i]f-'' Fig. 16.51 shows the MC history of the 
normalized [0i]f for the 16 replica of run [A2]. Here r refers to the molecular 
dynamics time for each replicum. While on this scale the bulk of the data are 
below one and hence not visible, there are several peaks, which give a significant 
contribution to the mean value. These spikes also affect the error estimate 
through both the variance and the integrated autocorrelation time. For statisti- 
cally accessible quantities the error should approach a 1/\/t behavior in the limit 
r^oo. In this respect we found /i and all other fermionic correlation functions 
we considered to be very hard to measure. Even when averaging over 16 replica, 

^Prom cq. fRTa|l it follows that [0i]f = Ti^St^t/ZL^. 




Figure 6.6: Normalized [(/)i]f and smallest eigenvalue from one "sick" replicum of 
run [A2]. Evidently the spike in [i;^i]f is dominating the statistical error cr(/i). 



this asymptotic behavior does not set in after r~500. 

The reason is the rare occurrence of very large values of [0i]f, which appear 
to be correlated with small eigenvalues of . However, this effect is washed out 
by using several replica. We therefore show in Fig. 16.61 the MC history of [0i]f, 
Amin and our error estimate for /i for one replicum of run [A2] with such a spike 
in [0i]f. Indeed, for each spike in [0i]f the smallest eigenvalue drops below its 
average. That the converse is not true could be ascribed to a lack of overlap 
of the eigenvector corresponding to Amin with the source needed to compute the 
quark propagator. Quantitatively, for the correlation between [(/)i]f and Amin 
we measure a value of C*[(^^]p^Amin = —0.33(4) if we use all replica and —0.46(6) 
from the replicum shown in Fig. 16.61 alone. Here we used as a definition of the 
correlation Ca,b between two (primary) observables A and B 



Ca,b = / \ /\ / ^ ^ ^Yiat - 1 < Ca,b < 1 . (6.16) 

Even though in the limit of infinite statistics configurations carrying very 
small eigenvalues are given the correct weight, depending on the algorithm this 
might be badly approximated for a typical ensemble size. Similar arguments 
referring in particular to the HMC algorithm motivated the introduction of the 
polynomial hybrid Monte Carlo (PHMC) algorithm in refs. [87,93,94]. 

Hence the difficulty in measuring fermionic correlation functions might be an 
efficiency problem related to the choice of the algorithm. To check this conjecture 
we employ a second algorithm and compare ensembles generated by HMC (with 
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two pseudo-fermion fields) with PHMC ensembles. Indeed, PHMC can be tuned 
in such a way that it enhances the occurrence of configurations carrying small 
eigenvalues, thus resulting in a better sampling of this region of configuration 
space. A reweighting step is introduced to render the algorithm exact. As a 
preparation for the following discussions we want to recall some properties and 
introduce the notations concerning the PHMC. 

6.6.1 The PHMC algorithm 

One of the main ideas of the PHMC algorithm is to deliberately move away from 
importance sampling by using an approximation to the fermionic part of the 
lattice QCD action. More precisely, in an HMC algorithm the inverse of is 
replaced by a polynomial Pn,e{Q'^) of degree n. Here Pn,t{x) approximates 1/x in 
the range e < x < 1. As a consequence this algorithm stochastically implements 
the weight VU det P~l{Q'^)e~^^ , whereas standard HMC generates ensembles ac- 
cording to VU det Q^e~^o . Denoting averages over the PHMC ensemble by (. . .)p, 
the correct sample average of an observable (O) can then be written as 

W 

{O) = {Ouj)p , where u = , (6.17) 

and we introduce the reweighting factor as a (partially) ^ stochastic estimate of 
det{Q^-Pn,e(Q^)}- When using Chebyshev polynomials the relative approximation 
error of P„^e is bounded by 5 ^ 2exp(— 2-^?^) in the range e < x < 1 . 

To give an impression of the role of e and 5 we plot in Fig. 16.71 a set of poly- 
nomials Pn,t{x) for typical (in our simulations) values of these parameters and 
compare them with l/x in the region of small x. Depending on the smallest 
eigenvalue of the parameters e and n have to be tuned such that the reweight- 
ing factor does not fluctuate too much. The authors of ref. [87] suggested to take 
e of the same order as (Amin) and in practice used e ~ 2(Amin) and 5 < 0.01. 

Recalhng that PHMC replaces detQ^ in the HMC weight with detP~^^((0^) 
and observing from Fig. 16. 71 that Pn,t{x) is smaller than 1/x for x < e, the afore- 
mentioned property of enhancing the occurrence of small eigenvalues is evident. 
At this point we would like to note that the fermionic contribution to the driv- 
ing force in the PHMC is bounded from above since Pn,e{x) is finite even at 
X = 0. In this way the polynomial provides a regularized inversion of Q^, thus 
also addressing the problems mentioned in Section 

6.6.2 HMC vs. PHMC 

Coming back to the comparison of samples from HMC and PHMC, we repeated 
run [A2] with PHMC using a polynomial of degree 140 and e = 6 ■10"'', resulting 

^Through the separate treatment of the lowest eigenvalues of Q'^ the infrared part of W is 
evaluated exactly. 
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Figure 6.7: Three different Chebysliev polynomials approximating 1/x, all with 
6 = 0.001. The right-hand plot shows the relative deviation from 1/x as a function 
of X. There the points (e, S) are marked by dots. 



in 5 ~ 0.002. The ratio e/ (Xmin) turned out to be around 2.7. In Fig. l6.8l we plot 
for this run the MC history of [4>i]f and of [</)i]f ■ ^, which enters into eq. ()6.17jl 
if we consider O = [4>i]f, i.e. 

/. = ([^,l,> = ([^j,..)p=<I:*^|±l£. (6.18) 

We first observe that apart from removing the largest spikes the inclusion of the 
reweighting factor does not seem to significantly change the relative fluctuations. 
This means that the parameters of the polynomial have been chosen properly. 
Events where [0i]f assumes a value 0(10) times larger than fi are no longer 
isolated as in Fig. Ifi.5l but happen frequently, which means that the PHMC algo- 
rithm can more easily explore the associated regions in conflguration space. This 
is what allows a reliable error estimate as shown in the upper part of Fig. 16. 8[ 
i.e. with 16 replica the asymptotic behavior of the error sets in after r^lOO. 

The advantage of using PHMC instead of HMC can be clearly seen by con- 
sidering the spread of (j{fi)\/T among the replica. We analyzed this quantity in 
extensions of runs [A2] and [A3] . The result is shown in flgure Fig. 16.91 where the 
shaded areas represent the range of values covered by the 16 replica as a func- 
tion of the MD time. In the limit of inflnite statistics all replica should converge 
to the same value, which need not be the same for the two algorithms because 
of reweighting and different autocorrelation times. We see that the spread for 
the HMC data is more than twice as large as for PHMC, i.e. the error on fi is 
signiflcantly harder to estimate with HMC. 

What we are suggesting here is that the algorithm should be chosen depending 
on the type of observables and the parameter values. From our experience we 



50 



Algorithmic issues 



0.04 |-r 




100 200 300 400 500 

r 



Figure 6.8: Monte Carlo history for the 16 replica of run [A3] showing the correla- 
tion function [0i]f and the product [i;^i]f -u;, where oj is the normalized reweighting 
factor. Our error estimate of /i shows the expected scaling behavior as soon as 
the run is long enough for a reliable extraction of Tint- 



conclude that PHMC sampling might just be more effective than HMC when 
computing fermionic quantities that are sensitive to small eigenvalues. 

To gain some more insight into the difference in sampling we consider the 
distribution of Amin since this is where we expect the largest effect. The distri- 
butions are analyzed by treating Abin = Xbin(Amin) as a primary observable. Here 
Xbin denotes the characteristic function of each given bin in the histogram. We 
then perform our normal error analysis for (Abin), where eq. fj6.17p has to be used 
if it is a PHMC sample. For comparison (Abin) p is also analyzed in this case. 

The histograms in the upper part of Fig. 16.101 compare the results from 200 
independent measurements produced by HMC and PHMC (runs [A2] and [A3], 
respectively). As expected the distributions agree within errors. For the PHMC 
run we also plot the unreweighted histogram, i.e. (Abin)p- Here we again con- 
firm that with the parameters we chose for the polynomial the PHMC produces 
more configurations with small eigenvalues than HMC. As a consequence of the 
reweighting the errors at the infrared end of the spectrum should be smaller for 
the PHMC data. This is explicitly verified in the lower part of the plot where we 
show the ratio of the errors on (Abin) from the two algorithms. One can see that 
below ~ 7-10~^ the error from PHMC is at least a factor two smaller than from 
HMC, with a more pronounced difference towards even smaller eigenvalues. The 
advantage in using PHMC to sample this part of the spectrum is significant and 
we will make use of this in the following discussion. 
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Figure 6.9: Monte Carlo history of a{fl)^/T from extensions of runs [A2] and [A3]. 
For the two algorithms we show the ranges covered by the 16 replica. 



6.7 Comparison to the quenched case 

In the previous section we studied various problems related to the occurrence of 
small eigenvalues. All the data presented there were produced at bare parame- 
ter values, which correspond to relatively large quark masses and small volumes. 
These small eigenvalues might therefore have a different nature from the "physi- 
cal" ones expected to show up in large volumes and/or close to the chiral limit. 
Here and in the next section we will establish them as cutoff-effects. 

To this end we made an additional simulation at the parameters of run [A2] 
and calculated the ten lowest-lying eigenvalues Aj, i = 1 ... 10. In Fig. 16.111 the 
smallest eigenvalue, Ai, is denoted by an open symbol. It seems that while A2 
through Aio form a rather compact band, the lowest eigenvalue fluctuates to very 
small values quite independently of the others. It is expected and has been shown 
numerically [95] that the spectrum of the Dirac operator depends quite strongly 
on the bare gauge coupling. A well-defined lower bound should be recovered close 
to the continuum limit only. Therefore we take the strong fluctuations of Amin as 
an indication for the presence of large cutoff-effects. Here we should point out 
that the eigenvalues of the Dirac operator are not on-shell quantities and hence 
the Symanzik improvement programme does not necessarily reduce cutoff-effects 
here. Quenched experience even suggests that the opposite might be true [96]. 

The occurrence of small eigenvalues at these bare parameters poses a some- 
what unexpected problem in dynamical simulations. Comparing the quenched 
situation to the Nf = 2 dynamical case, the naive expectation is that at fixed 
bare parameters the probability of finding configurations with small eigenvalues 
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Figure 6.10: Upper panel: histograms of Aminj i-e. (Abin) vs. 'bin', from runs [A2] 
and [A3]. For run [A3] we also show (Abin)p- Lower panel (from the same data): 
ratio of the error on (Abin) from HMC to that from PHMC. 



should be reduced by the determinant. To us the more relevant question seems 
to be whether small eigenvalues are suppressed in a situation where the physical 
parameters (e.g. volume and pseudo-scalar mass) are kept constant. 

Using the quenched data from ref. [72] and the dynamical data from refs. [97] 
and [98] (where an estimate of ro/a = 5.21(6) for /3 = 5.2 can be found) we chose 
the parameters of the quenched run [A5] such that the lattice spacing and the 
(large volume) pseudo-scalar mass are matched to run [A4]. This was found to 
occur at almost equal bare current quark mass (see Lm in Table IFTTj) . In Fig. l(i.l2l 
we compare the distributions of Amin for these two runs. Two comments are in 
order here: 

• For the dynamical run the mean value is shifted up from 1.44(1) -10"^ to 
1.72(5) ■ 10~^. This agrees with the naive expectation but in a physically 
matched comparison it is a non-trivial observation. 

• The distribution itself is significantly broader compared to the quenched 
case and in particular it is falling off more slowly towards zero. This means 
that even though (Amin) is larger for Nf = 2 the probability of finding very 
small eigenvalues is enhanced. 

The second point, i.e. that the lower bound of Amin is less well-defined, seems 
to imply that at a lattice spacing of a ~ 0.1 fm the cutoff-effects are much larger 
in the Nf = 2 case. To substantiate this we will compare the distribution of Amin 
to that from a run at finer lattice spacing and matched physical parameters. 



6.8 Finer lattices 



53 



A 



10 



-3 








100 



200 



300 



400 



measurement 



Figure 6.11: Monte Carlo history of the ten lowest eigenvalues at the parameters 
of run [A2]. The open symbols denote Amin- 



6.8 Finer lattices 

Apart from cutoff-effects, in the massless theory the Schrodinger functional cou- 
pling g"^ is a function of the box size L only [9, 10]. We measured it on a small 
lattice of extension L/a = A at /? = 5.2, obtaining a value of ^^ = 3.7(1). We then 
extrapolated to this value the L/a = Q data used in ref. [99] as a function of {3. Our 
result from the matching is that for the two-flavor theory a bare gauge coupling 
oi (3 = 5.5 roughly corresponds to a lattice spacing, which is 1.5 times smaller 
than at /3 = 5.2. This estimate will be confirmed by other non-perturbative data 
as well as the perturbative evolution of a in Chapter [7| 

Hoping that the algorithmic difficulties arising from cutoff-effects would be 
much smaller in this situation, we simulated a 12^ x 27 lattice at this value of 
/3 (run [A6]) using the HMC algorithm. With the k we chose (and ignoring the 
change in renormalization factors) the bare quark mass Lm is roughly matched 
to the heavier runs at /5 = 5.2. We therefore compare run [A6] with run [A3]. 

Normally, a constant acceptance requires a decrease of the MD integration 
step-size if ones goes to finer lattices at fixed physical conditions. This argument 
is based on the scaling of the small eigenvalues^, which influence the MD driving 
force. We found that (Amin) in run [A6] is a factor two smaller than in run [A3]. 
Nevertheless, at /5 = 5.5 the step-size necessary for a certain (^ 90%) acceptance 
is roughly the same as at /3 = 5.2. This indicates that the value of 6t we had to 
use in the HMC runs aX (3 = 5.2 was dictated by the occurrence of extremely small 
eigenvalues rather than by the average smallest eigenvalue. In addition, where 

^As the squared Dirac operator, its eigenvalues have mass dimension two and should thus 
behave like a^. 
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Figure 6.12: The histograms of Amin from run [A4] {Nf = 2) and run [A5] {Nf=0). 
Despite its higher mean value the dynamical data show a much larger probability 
of finding very small eigenvalues. 



in run [Al] at the same average acceptance a maximum of 47 proposals were 
rejected in a row, the maximum for run [A6] is 4 trajectories. For this reason 
e~^^ shows no autocorrelation and its distribution is well separated from zero. 

Concerning fermionic observables, we have not observed spikes and hence 
expect the error to scale properly. However, for an accurate estimate of the error 
on e.g. /i our present statistics is not sufficient. One should note that ratios of 
correlators are easier to estimate since usually numerator and denominator are 
correlated, which reduces the impact of statistical fluctuations. This applies to 
essentially all relevant quantities, including the axial current improvement and 
renormalization constants discussed in the next chapters. 

The reason for these effects is the change in the distribution of Amin- To 
compensate for the different lattice spacing. Fig. Ifi.lHl compares Amin/(Amin) from 
runs [A3] and [A6]. One can clearly see that at the flner lattice spacing the 
probability of flnding a smallest eigenvalue less than half its average is greatly 
reduced compared to (3 = 5.2. The width of the distribution is smaller in this case 
and in particular the spectrum is now clearly separated from zero. Quantitatively, 
the normalized variance of Amin is reduced from 0.18(1) to 0.13(2). 

This comparison explicitly shows that the long tail of the eigenvalue distribu- 
tion we observed at a ^ 0.1 fm, and which caused the problems we have discussed, 
is a cutoff-effect. Matching also run [A6] to a quenched simulation (run [A7]), 
we again found an upward shift of (Amin) for the dynamical case. In addition, at 
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Figure 6.13: Normalized distributions of Amin from runs [A3] (/3 = 5.2) and [A6] 
(/3 = 5.5). While the data from the coarse lattice stretch almost to zero, the (3 = 5.5 
data seem to have a more well-defined lower bound. 




Figure 6.14: The histograms of Amin from run [A6] {Nf = 2) and run [A 7] (iVj =0). 
At this finer lattice spacing the lower end of the spectrum appears to be similar in 
the quenched and the dynamical case. 
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this finer lattice spacing, the tails of the distributions of Amin look already very 
similar to each other as shown in Fig. 16. 141 

6.9 Conclusions 

At a lattice cutoff of approximately 2 GeV, corresponding to a^O.l fm, we have 
studied the behavior and performance of HMC-type algorithms in an interme- 
diate size volume of Ifm^. We discussed problems related to the occurrence 
of small eigenvalues in two-fiavor dynamical simulations with improved Wilson 
fermions. We found these small eigenvalues to be responsible for large Hamilto- 
nian violations in the molecular dynamics. Even for integration step-sizes such 
that the acceptance is 80 ~ 90% those can still cause long periods of rejection, 
thus degrading algorithmic performance. However, in spite of employing only 
single-precision arithmetics we never observed reversibility violations. 

In addition, those eigenvalues make the estimate of fermionic quantities very 
difficult. The naive intuition is that the fermionic determinant should suppress 
small eigenvalues compared to the quenched case. Through a direct comparison 
at matched physical parameters we indeed verified that (Amin) is larger with two 
dynamical fiavors. On the other hand there is no obvious expectation for the 
tail of the distribution and we observed that it extends further towards zero 
than in the quenched case. Given the infrared cutoff induced by the Schrodinger 
functional boundary conditions and the quark mass we interpret this as a lattice 
artifact. We were able to confirm this picture with a simulation at finer lattice 
spacing, where the spectrum turned out to have a much sharper lower bound. 

In our study we found that the PHMC algorithm is more efficient than HMC 
(with two pseudo-fermions) in incorporating the contribution to the path integral 
of configurations carrying small eigenvalues. In other words, the distortion of 
the spectrum by cutoff-effects actually makes it advantageous to deviate from 
importance sampling. Also without such special problems we found PHMC at 
least comparable in performance to HMC (in our implementations). 

We want to emphasize that the problems discussed here do not occur only 
in the Schrodinger functional setup. Without this infrared regulator they are 
expected to show up already at larger quark masses. 

Due to the in a sense "more robust" nature of PHMC one might ask why 
one should not use it in all simulations. The reason is that with decreasing 
lattice spacing the smallest eigenvalues will also decrease, which in turn requires 
a smaller e in the construction of the polynomial as otherwise the fluctuations of 
the reweighting factor will render the algorithm inefficient. 

However, a smaller lower bound in the polynomial approximation implies a 
larger degree of the polynomial (at constant relative deviation 6). In the current 
implementation of the PHMC, single precision restricts the degree to ~ 140 since 
otherwise roundoff problems appear in the construction of the polynomial. 
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We therefore used the PHMC only in the simulation for the axial current 
renormalization at j3 = 5.2 and 5.29. The simulations for the axial current im- 
provement condition are performed with a smaller time extension and larger quark 
mass. Therefore one can safely use HMC at all lattice spacings we consider. 



Chapter 7 

Axial current improvement 



In the quenched approximation, the improvement coefficients Cgw and ca have 
been determined non-perturbatively in the relevant range of bare couphng (or 
lattice spacings) in [12]. The improvement conditions, from which they have 
been determined, were derived from the chiral symmetry of the theory in the 
continuum limit. 

More precisely, the PCAC relation ()2.28j] was required to hold at finite lattice 
spacing [12,61,63]. As will be detailed in the next section, the PCAC relation is 
considered considered with different external states. In [12,61,63], finite volume 
states were chosen, formulated in the framework of the Schrodinger functional. 
Later, a determination of ca was performed by evaluating the PCAC relation in 
large volume [100,101] at several values of the lattice spacing. While at a ^ 0.1 fm 
these results for ca differ quite significantly from the finite volume definition [12], 
at smaller lattice spacings the difference decreases. 

For the interpretation of this difference, one should keep in mind that be- 
yond perturbation theory the improvement coefficients themselves are affected 
by 0(a) ambiguities. In some detail this has been discussed and demonstrated 
numerically in [102]. The 0(a) ambiguity simply corresponds to the fact that 
in the implementation of the improvement programme the theory is treated only 
up to O(a^) effects. Thus any change of 0(a) of the improvement coefficients 
only contributes to the terms 02 and £2 in fl3.14| I3.15j) and does therefore not 
invalidate 0(a) improvement. 

While this obviously forbids a unique definition of the improved theory, the 
0(a) ambiguities can be made to disappear smoothly if the improvement con- 
dition is evaluated with all physical scales kept fixed, e.g. in units of tq fl3.13|) . 
while only the lattice spacing is varied [102]. This is what is meant by a "line of 
constant physics" (LCP). 

Since Symanzik's effective theory describes the lattice artifacts only asymp- 
totically (Section 13. 3^ it can only be valid for matrix elements dominated by 
states with energy E ^ a~^. It is therefore important to make sure that the 
improvement conditions are also imposed using low energy states. So far, the 
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methods of [100, 101] have not yet been implemented such as to satisfy these 
conditions. 

In [103], two improvement conditions for ca, were studied, for which one can 
choose the kinematic parameters such that the above criteria are satisfied. They 
are formulated in finite volume with Schrodinger functional boundary conditions, 
which furthermore helps to render the numerical evaluation feasible in full QCD. 
These improvement conditions will be discussed in the next section, where we 
motivate our choice of one of them to compute ca in the Ni = 2 theory. 



7.1 Strategy and techniques 

Before going into the details of our strategy and techniques let us comment again 
on the constant physics condition. In order to keep the physical volume constant, 
we need to know how the lattice spacing depends on (3. With this knowledge we 
can tune (3 such that a certain L/a corresponds to a prescribed value of L (or 
L/tq). However, this tuning of the volume is not critical in the evaluation of 
Ca, since the latter depends on the volume only through effects of order a/L} A 
relative uncertainty A in the physical volume L (or equivalently in ro) thus implies 
a relative uncertainty in ca, which is proportional to ajL-!S.. As a consequence, 
even if A varies a bit in the considered range of lattice spacings, this is quite 
irrelevant, in particular if A is a smooth function of the lattice spacing. Therefore 
the constant physics condition with respect to the volume has to be enforced with 
only moderate precision. 

In the remainder of this section we will state in more detail how the constant 
physics condition is implemented. We will also discuss the methods in [12, 103] 
to determine ca, with emphasis on the one we finally used. 



7.1.1 Constant physics condition 

With two degenerate flavors of light quarks, the theory has two bare parame- 
ters, (5 and m^. Roughly speaking, the bare quark mass mg controls the physical 
quark mass and the bare coupling determines the lattice spacing, defined at van- 
ishing quark mass (for a more precise discussion see [63,85]). Non-perturbative 
estimates of 

[ro/a](5.2) 
[ro/a](/5) 

are available in a limited range of /3 [97,104]. In [2], the results of [97,104] were 
extrapolated to zero quark mass. Taking directly these values for [ro/a](/3), we 
have the points with error bars in Fig. 17.11 From those we roughly estimated the 
location of the filled points, using the perturbative dependence of the lattice spac- 
ing on /? as a guideline. For our action this is given to three loops in [105], which 

^Effects of order am will be discussed later. 
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Figure 7.1: The evolution of the lattice spacing a with the inverse bare gauge 
coupling (3 from perturbation theory in the lattice scheme and large volume data 
[104] for the scale tq. The filled points correspond to our "scaled" simulations. 

builds on various steps carried out in [106-112]. Applied as a pure expansion in 
the bare coupling (no tadpole improvement), one has 



The evolution of the lattice spacing relative to our reference point at {g'^)'^ = Q/5.2 
is then expressed by the function t(/5) = a(6//3)/a(6/5.2), which is plotted as a 
thick line in the graph. It confirms that the filled points are very reasonable 
choices. Note that other forms of applying bare perturbation theory (differing 
from eq. ()7.2|) in the (y'Q-term) would give somewhat different results, but since we 
are interested in a rather limited range in gQ this is of no great importance. Later 
we will show that systematic uncertainties in ca introduced by this approximate 
scale setting are negligible. Finally, we keep the PCAC mass ()5.22j) . evaluated 
with an average of fx and gx approximately constant. 

7.1.2 Improvement conditions for the axial current 

In this section we discuss criteria for the choice of the improvement condition. 
Quark masses derived from the PCAC relation are of the form 



a(^o) 



e-k7^-(.i)-]/2^o[^2/(^,)2]-W2.g [l + qlgl- (,^)2] + O {{g',y) ] , 



aim 




■m{x] a, (3) 



2{a\P^{x)\(3) 



(7.3) 
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Since tliis mass is obtained from an operator identity, it is independent of the 
states |a) and \I3) as well as the insertion point x up to cutoff-effects. Enforcing 
this independence at finite lattice spacing leads to possible definitions of improve- 
ment conditions [103]. Inserting the expression for the improved current ()4.25p 
in the previous equation, the quark mass can be written as m = r + ac^s + O(a^) 
with 



If we now consider two sets of external states and two insertion points, the im- 
provement condition m{x\ a, (3) = m{y; 7, 6) yields 



The "sensitivity" to ca of an improvement condition is therefore given by aAs, 
since it is this term, which is multiplied by ca in the current quark masses. 

Once a reasonably large sensitivity is achieved, all improvement conditions 
at constant physics are equally valid in the sense that 0(a) effects are removed 
in on-shell quantities. However, the way in which higher-order lattice artifacts 
are modified will depend on the concrete choice of the improvement condition. 
In particular, if states with energy not so far from the cutoff are involved, large 
O(a^) effects might be introduced. 

To discuss the different improvement conditions we need to generalize the 
Schrodinger function correlators introduced in Section EUl Instead of projecting 
both boundary quarks to zero momentum separately ()5.13l5.14p . we now use 
boundary fields products, which depend on spatial trial "wave functions" uj{x.) 
[113]. In the following we will use 



and 




(7.4) 



(7.5) 



- Ca = 



aAs a s{x]a, (3) — s{y;j,6) 



(7.6) 




(7.7) 



X 




(7.8) 



X 



and 




(7.9) 



which depend on the pseudo-scalar operator 




(7.10) 
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at the Xo = boundary and the corresponding operator 0''^{uj') at the upper 
boundary of the SF cyhnder. 

The Schrodinger functional version of eqs. ()7.4|) and (|7.5p is then given by 



To determine ca in the A'f = theory [12], Ar and As were originally defined 
through a variation of the periodicity angle 6 of the fermion fields, while keeping 
xo = T/2 and uj = const fixed. For this method the sensitivity aAs is quite low 
when L > 0.8 fm, T = 2L. In addition, with dynamical fermions different values 
of 6 would require separate simulations. We therefore consider this method as 
too expensive and disregard it in the following. In the quenched approximation 
two alternatives have been explored in [103]. 

Requiring the quark mass to be independent of xq (for fixed 6 and u = const) 
is technically easy to implement. However, also in this case the sensitivity is 
small unless large values of 6 are used. Moreover, the contribution of excited 
states is not well controlled, because one insertion point must be rather close to 
a boundary to achieve a sufficiently large sensitivity. Thus, energies which are 
not far removed from the cutoff may contribute. 

The simultaneous requirements of large sensitivity and control of excited- 
state contribution can be fulfilled more easily with the second method, where 
variations of the wave function u have been considered. Ideally, one would like 
to use two wave functions u^.^ and uJt^^i, such that the corresponding operator 
0"'{u!) couples only to the ground and first excited state in the pseudo-scalar 
channel, respectively. As one can easily see from eq. ()7.12|) the sensitivity to ca is 
then proportional to ^ —fnl c where m,r,n denotes the mass of the nth excited 
state in the pseudo-scalar channel. Higher excited states are (by definition) not 
contributing and in principle the method can be used for rather small T. Hence, 
we find this the most attractive method both from a theoretical and practical 
point of view. In the next section we will detail our approximation to this ideal 
situation. 



7.1.3 Wave functions 

We will now proceed to the more technical aspects of our method. In order to 
approximate cu^^o and cj^^i consider a set of N wave functions. Given a vector u in 
this A-dimensional space, projected correlation functions are defined as {u, /a) 
and {u,fiu), i.e. /a is regarded as a vector and fi as a matrix in this space. It 



and 






(7.12) 



2fp{xo]u) 
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is useful to represent fx (X = A, P) and /i as [114] 

M-l 

/x(xo; uj,) = J2 F>f e-™--"^o + 0(6-"^-'"^°) + 0(e-"^G(T-xo)) ^ ^^ ^^^ 

n=Q 
M-l 

h{uj,,ujj) = ^ i;>;e-™--"^ + 0(6-"^-^^^) + 0(e-"«^) , (7.14) 

n=0 

where n labels the states in the pseudo-scalar channel in increasing energy and 
is the overlap of such a state with the one generated by the action of C"(a;j) 
on the SF boundary state. The mass mc belongs to the lowest excitation in the 
scalar channel, the O"*"^ glueball and the coefficients are proportional to the 
decay constant of the nth state. Here we have suppressed the explicit volume 
dependence of all quantities. 

Knowledge of would allow the construction of vectors m", such that - up to 
corrections of order e~"^^''^''^ - the correlation (m", /a) receives contribution from 
the nth state only. These may be computed from the f " by a Gram-Schmidt 
orthonormalization. Clearly, m° and can then be used to approximate cj^^o and 

An approximation to the can be obtained from the eigenvectors of the 
positive symmetric matrix /i. We denote the normalized eigenvectors of /i by 
?7°,?7\ . . . corresponding to eigenvalues A° > > . . . and apply eq. ()7.14jl to r]^ 
to obtain^ 

v°(^;°, ryO)e-™-^«^ + 0(e-'"-i'^) = AV . (7.15) 

Clearly Aq is of order |woPe~'"'"^''^, which is inserted to rewrite the previous equa- 
tion in terms of the normalized vectors = v^/\v'^\ 

yO({;0, r/°) = + 0(e-('"-'^-™--«)^) . (7.16) 

Since {)° and rj^ are normalized vectors, it follows that {v^,ri'^) = 1 up the error 
term given above. Together with the orthogonality of the eigenvalues of /i this 
implies 

11^)0 -r/0||2 = 0(e-("^-^-'"-°)^) (7.17) 
and (r/\{)°) = 0(e-(™--^-'"-°)^) . (7.18) 

Thus, to the order indicated above, -O^ is given by rj'^ and rj^ is orthogonal to 
the "ground state vector" -0°. As eigenvectors of a symmetric matrix the r/" are 
already orthogonal and we therefore use the approximation 

^^fi-'^Vi^i and cj^^i ~ ^ r/^^cjj (7.19) 

i i 

^The glueball contribution will be dropped from now on. 
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to obtain correlators, which are (for intermediate xq) dominated by the ground 
and first excited state, respectively. We note in passing that the ratios v^/v^ 
have a continuum limit if the wave functions are properly scaled with the lattice 
spacing. 

In our simulations we restrict ourselves to a basis consisting of three (spatially 
periodic) wave functions defined by 

t^i(x) = A^"^ ^^i(|x-nL|), i = 1,...,3, 

neZ3 

uj.,[r) = r-^l\-r/{-i-o) ^ (7,20) 

where is some physical length scale. We thus keep it fixed in units of L, choos- 
ing qq = L/6. The sum over n is required to preserve the spatial periodicity. 
In practice the summation is stopped when the norm of u no longer changes 
within single precision arithmetics. The (dimensionless) coefficients Ni are fixed 
to normalize the wave function via '^^^fi'^) = 1- In this context the correla- 
tion functions introduced in Section 15.21 can be regarded as belonging to the fiat 
wave function u;o(x) = where both quarks are projected to zero momen- 

tum separately. In this case x and y in eq. ()7.10j) are uncorrelated and thus full 
translational invariance can be used without performing additional inversions of 
the Dirac operator (see Appendix IB.l.lj) . In the general case uJi,,,^ we replace 
one of the spatial sums in eq. ()7.10|) by a sum over eight far separated points, 
which means that one performs eight times as many inversions. This additional 
computational effort is still small compared to the cost of the inversions in the 
HMC update. 

7.2 Numerical computation 

7.2.1 Results for the improvement coefficient 

Table 17.11 summarizes the parameters of our simulations for the axial current 
improvement constant. The values of the other improvement coefficients are 
given in Sections 14.21 and 15.11 

The (3 values for run [C2] and [C3] have been chosen such that L/tq is approx- 
imately the same as in run [CI], which corresponds to L ~ 1.2 fm. In exploratory 
quenched studies [103] this volume was found to be sufficient for the described 
projection method to work. N^cas is the number of estimates of ca (with the num- 
ber of replica denoted explicitly), separated by Tmeas in Monte Carlo time. The 
autocorrelation of these measurements turned out to be negligible. The column 
labeled am/t{[3) refers to the bare quark mass m = r{T /2]ujQ)+acxs{T /2]ujq), cf. 
eqs. ()7.1H l7rT2|) . which is equivalent to ()5.22j) . The 1-loop value of ca from [61] 
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12 


12 
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5 
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16 
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8 
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0.0366(36) 


[C6] 


24 


24 


6.12 


0.136139 


2x21 


12 


0.0002(4) 


0.0244(21) 



Table 7.1: Summary of simulation parameters and results for ca- Runs [Cl]- 
[C3] are at constant physics. 



is used here. We tuned the hopping parameter k in order to keep am/t{l3) fixed 
when varying (3, thus ignoring (presumably small) changes in the renormalization 
factors. Note that we have chosen a finite, but small bare quark mass of around 
30 MeV. Such a mass helps (in addition to the Dirichlet boundary conditions) to 
reduce the cost of the simulations. Results from the remaining simulations are 
used to discuss systematic uncertainties in our determination of ca- 

We employed the hybrid Monte Carlo algorithm with two pseudo-fermion 
fields as described in Section lUTTl For all observables we have checked the expected 
scaling of the statistical error with the sample size and thus verified the absence 
of the problems described in Sections 16.51 and 16.61 at the volumes and masses we 
consider here. For run [C2] Fig. 17.21 shows how the error of our estimate for 
Ca (see below) behaves as a function of the Monte Carlo time r. One can see 
that already after approximately 200 trajectories the region of statistical scaling 
is reached. 
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Figure 7.2: The error on c\ vs. the Monte Carlo time for run [C2]. 
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Figure 7.3: The effective mass in lattice Figure 7.4: As(xo) and Ar(xo) deter- 
units of tlie projected correlation func- mined from lOt^^ and ujjr^i in run [C2]. 
tions {r]^, fp) and (r/^, /p) from run [C2]. 



In Fig. 17. HI we show the effective masses from /p(xo; cJtt.q) and /p(xo; cuTr.i) as 
obtained in run [C2]. Two distinct signals are clearly visible, which indicates that 
the described approximate projection method works well at these parameters. 
The energy of the first excited state is not far away from a~^, suggesting that in 
even smaller volumes the residual O(a^) effects would grow rapidly. In the spirit 
of the remark after eq. ()7.19j) at the other values of /5 we used the same linear 
combination of wave functions to define Ut^^q and Wj^^i, namely 

77° = [0.5172, 0.6023, 0.60811 
' 7.21 
and r/i = [0.8545,-0.3233,-0.4066], 

which are the ones determined in run [C2]. When scaled in units of ro, this yields 
effective masses similar to those shown in Fig. l7.3l Results from a redetermination 
of 77^") and T]^^^ in the other matched simulations agree with eq. ()7.2H) to a high 
precision. In fact, one can easily distinguish the matched and unmatched runs 
using e.g. 77°. Comparing the two 24^ runs 

[C3] 7/° = [0.5173(2), 0.6024(1), 0.6079(1)] , (7.22) 
and [C6] r/° = [ 0.5126(2), 0.6042(1), 0.6101(1) ] , (7.23) 

it is clear that the eigenvector rj^ from run [C3] is in good agreement with ()7.21|) . 
while a significant deviation is seen in run [C6]. 
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In Fig. 17.31 the error on the effective mass of the first excited state is seen to 
be quite large, but what actually enters the computation of ca is the error of 

Ar(xo) = r{xo;uj^^i) -r{xo;uj^^o) (7.24) 
and As(xo) = s{xq;uJt,^i) — s(xo;u;^,o) • (7.25) 

These profit from statistical correlations of the correlation functions entering their 
definition and thus have smaller statistical errors as can be seen in Fig. l7.4t where 
we plot aAr and a^As from the same data used in Fig. 17.31 

Fig. 17.51 collects results for the "effective" ca{xo) = — Ar(xo)/aAs(xo) from 
the matched runs [C1]-[C3]. We see little variation for xq ^ 6a, which we take 
as another signal that high energy states, which could contribute large 0(a) am- 
biguities in the improvement condition, are reasonably suppressed in this region. 
We complete our definition of ca with the choice Xq = T/2, which is at the same 
time scaled in physical units and in agreement with the xq ^ 6a bound for all 
our lattices. 

Finally, ca is plotted as a function of Qq in Fig. 17.61 The solid line is a smooth 
interpolation of the data from the matched simulations, constrained in addition 
by 1-loop perturbation theory: 

. o 1 - 0.4485 

c^toS = -0.00756 ,,2 x^-^^. (7.26) 

It is our final result, valid in the range 0.98 < Qq < I-IQ within the errors of the 
data points (at most 0.004). 
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Figure 7.5: Effective c\ as determined from lo-j^^q and cj^r,! for runs [C1]-[C3]. 
Points with xo/a = 6 are marked by filled symbols. 
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This non-perturbative result is quite far away form 1-loop perturbation the- 
ory, in particular at the coarsest lattice spacing. For (3 = 5.2, the perturbative 
value of ca = —0.0087 is almost an order of magnitude smaller. 

Of course this has a significant impact if simulation data are now analyzed 
using the non-perturbative ca instead of its 1-loop value. Using e.g. data from 
[115], we see that the effect on the result for the pseudo-scalar decay constant at 
this lattice spacing is as large as 10%. 



7.2.2 Uncertainties due to deviations from the 
"constant physics" condition 

We should check whether the volumes in our runs [C1]-[C3] are scaled sufficiently 
precisely or if systematic errors need to be added to the statistical ones on ca to 
cover possible violations of the constant physics condition. Table 1 shows that the 
bare PCAC mass has been kept constant to within about 10%. A renormalized 
quark mass would differ by the multiplication with a Z-factor, which is a slowly 
(namely logarithmically) varying function of a. Considering the restricted range 
of lattice spacings we cover, such factors can be safely ignored. Run [C4] is done 
with a quark mass which is more than twice as large as the one in run [CI], with 
otherwise identical parameters. The barely significant difference in ca confirms 
that the small deviations from the "constant mass" condition can be neglected. 

We also need to examine the impact of the uncertainty in L due to our per- 
turbative (or asymptotic) scaling of the lattice spacing. It has been argued in 
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Section 17.11 that the difference to a proper non-perturbative scaling is rather 
small. Also, estimating a possible change by comparing 3-loop to 2-loop and 
non-perturbative scaling, gives a deviation in t{l3) which is smaller than 10% in 
the whole range of Fig. 17.11 and thus the same maximum deviation applies to 
L/a. That this is negligible can be seen from a comparison of run [C5] to the 
interpolating formula ()7.26|) . In run [C5], L/a is 20% lower than the proper value, 
but ca does not differ significantly from the fit curve. As will be discussed in the 
next section, comparison with significantly smaller volumes at smaller (3 values 
can lead to enhanced O(a^) effects. 

Finally, by run [C6] we verify that the dependence of ca on the kinematic 
parameters disappears quickly when going to even larger values of (3. In this run 
we used gauge configurations from the calculation of Z-p [78]. Although those 
were produced at m = 0, ^ = 0.5 and a much smaller volume^, the resulting ca is 
only approximately two standard deviations away from our fit. 

7.3 Discussion 

For the 0(a)-improved action with non-perturbative Csw [59], we have determined 
the improvement coefficient ca for (3 > 5.2, which roughly corresponds to a < 
0.1 fm. The improvement condition was evaluated on a line of constant physics, 
which is necessary in the situation when 0(a) ambiguities in the improvement 
coefficients are not negligible. 

That this is indeed the case here can be seen from an additional run at /5 = 5.42. 
Decreasing L/a from 16 to 12 at constant T/a gives ca = —0.0559(21), which 
is about 25% larger in magnitude than its value on our line of constant physics. 
This can be understood from the fact that the energy levels in the pseudo-scalar 
channel increase if the volume is decreased. But in order to safely exclude large 
0(a) ambiguities, improvement conditions should only involve states with energy 
E <^ a~^. On this requirement we had to compromise more than we would have 
liked to do already in our LCP volume, where the maximum values for Ea are 
about 0.7. With larger energies in smaller volume it is only to be expected that 
the value of ca changes in trying to compensate the larger O(a^) effects. 

Although one could have improved the situation, i.e. made Ea smaller, by 
going to somewhat larger values of L (and T), this would have made the numerical 
computation much more expensive. 

Large O(a^) effects have indeed been found in the A'f = 2, 0(a)-improved 
theory [85] at (3 = 5.2 and these may well be related to the not so small 0(a) 
ambiguity in ca that we just mentioned. 

This can only be investigated further by studying the scaling violations in 
quantities such as F^^rQ after improvement and renormalization. The next step af- 
ter successfully implementing non-perturbative improvement for the axial current 

•^If one would extend i(/3) to this value, this would indicate a required lattice volume of 43^. 
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is therefore its (also non-perturbative) renormalization, which will be discussed 
in the next chapter. 

Clearly, the method described here may also be useful to compute ca in the 
three flavor case, where Cs„ is known non-perturbatively with plaquette and 
Iwasaki gauge actions [60,116,117]. 



Chapter 8 

Axial current renormalization 



After the successful non-perturbative improvement of the axial current, we can 
now turn to the problem of its normalization. In the quenched case [11], a 
non-perturbative renormalization of the isovector axial current was achieved by 
requiring certain continuum chiral Ward identities to hold at finite lattice spacing. 
These relate matrix elements of the axial and vector currents and thus link the 
normalization of the two. The local vector current is normalized by enforcing 
isospin symmetry, again in the form of an integrated Ward identity. 

In fact, the situation is similar to the calculation of ca- Also here the results 
(i.e. Za and Zy) are ambiguous due to cutoff effects and in order to obtain 
them as smooth functions of we need to evaluate our normalization conditions 
on a line of constant physics, keeping all length scales fixed. As discussed in 
Section \4.'A\ the normalization conditions have to be set up at vanishing quark 
mass since we want to implement a mass-independent renormalization scheme. 

While this is in principle possible in the Schrodinger functional, we have seen 
in Chapter IHl that cutoff effects make simulations at small quark masses and 
coarse lattice spacings difficult. Although this problem is addressed efficiently 
through the use of the PHMC algorithm, it is still desirable to formulate the 
normalization conditions such that they have only very little dependence on the 
quark mass. This is the point where we improve the method employed in the 
quenched case [11] by deriving a normalization condition from the axial Ward 
identity including the mass term. This new normalization condition naturally 
reduces to the previously used one in the chiral limit, but numerical results can 
now be more easily extrapolated in the quark mass. 

After deriving the normalization conditions for the isovector currents, the 
advantage of using the new condition is demonstrated through its chiral extrap- 
olation in a quenched example. Then the results for two dynamical flavors are 
presented and summarized in interpolating formulae. We also discuss systematic 
errors due to deviations from the line of constant physics and other variations of 
the normalization condition. 
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8.1 Continuum Ward identities 

In Section 12 .HI we considered Ward identities derived from the flavor chiral sym- 
metry of the continuum QCD Lagrangian. Here we will generalize these to matrix 
elements, where (part of) the operator insertion has support in the region, where 
the fleld transformations are performed. As a result eqs. ()2.27|) and ()2.28p will 
be modifled by terms containing the variation of this "internal" operator. 



8.1.1 VWI 

Again we consider only the case of two degenerate quarks flavors and pick a region 
TZ with a smooth boundary dTZ. Suppose Cint and Oext are polynomials in the 
basic flelds localized in the interior and exterior of this region, respectively. If we 
perform a fleld transformation with support in TZ, the variation of the external 
operator Ocxt vanishes and for the case of an isovector vector transformation it 
follows that 



where ^vC'int = '^"'{v) ^v^int if V denotes the point, where O lives. We now take 
10°' to be a constant in TZ. Using Gauss' law the integrated vector current Ward 
identity (VWI) is obtained 



((5eant)axt) = -/ da^(x)(F;(x)antaxt) . (s.i) 

' Jan ^ ' 

8.1.2 AWI 

The same construction is also valid for the axial current variation except that the 
mass term in ()2.19p cannot be written as a surface integral^ 



(5;Oint)axt) ™ - / dcr^(x) (A^(x)a„ta 
' Jan ^ 



2m / d4a;(P'^(x)a,,axt) • (8.2 



int*^cxt 

n 



The axial Ward identity (AWI, 18. 2|) will later be used for massless quarks with 
Oint equal to the axial current component A'l at some point yElZ. for this case 



^Due to additive mass renormalization, it is the current quark mass, e.g. (j5.22|l . rather than 
the bare quark mass, which appears in the integrated Ward identity. 
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we have 



Jan ^ 



ITT71 



abc 



(8.3) 



8.1.3 Euclidean proof of the Goldstone theorem 

As an immediate application of tlie integrated Ward identities witli operator 
insertions we will now show that a non-zero chiral condensate (ipip) leads to 
long-ranged current-density correlation functions [34], which in turn implies a 
massless propagating particle. This is the basis for the conjectured mechanism 
of spontaneous chiral symmetry breaking. 

We know from ()2.28|1 that for all x 7^ and zero quark mass the correla- 
tion function (c}^74^(a;)P"(0)) vanishes. From Lorentz invariance one can then 
conclude that the integrated version is of the form 

(A^(x)P'^(O)) = . (8.4) 

The constant k can be fixed by applying ()8.2j) with a sphere of radius r around 
the origin as the region 71, m = 0, Oi^t = P°'{0) and Ocxt = 1, which gives 

/ da,{x){A;ix)P-{0)) = -{6^P'^m = -K^^) . (8.5) 

J\x\=r 

In the last step the definition ()2.14j) has been applied to the pseudo-scalar density 
()2.20|) . Since the integrand's divergence is zero everywhere except at the origin, 
in an application of Gauss's law the entire contribution would come from the 
contact term at the origin. However, we can simply plug ()8.4j) into ()8.5|) and 
solve the trivial integral directly. The result is the surface of the 4-sphere times 
k such that 

- 1(7/;^) = k27r^ (8.6) 

and we finally arrive at 

(A^(x)nO)) = -A(^^)3i_. (8.7) 

Thus, if the chiral condensate (ipip) is non-zero, the correlation function on the 
left-hand side of ()8.7|) has no exponential decay and therefore the energy spec- 
trum of the theory has no gap. The massless particles propagating in the corre- 
lation function ()8.7|1 are the Goldstone bosons of the broken chiral symmetry. 
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8.2 Normalization conditions 

The normalization conditions for the isovector currents are obtained by deriving 
an identity between matrix elements in the continuum from the isovector sym- 
metries of the QCD action. The renormalized improved currents on the lattice 
are then required to satisfy these identities at finite lattice spacing. This results 
in conditions for the renormalization factors Zp^ and Zy, cf. eqs. ()4.34l H7^ . 

In a massless renormalization scheme these conditions have to be set up a 
vanishing quark mass. If numerical simulations at the critical point are not 
possible, the normalization factors become the result of a chiral extrapolation. 

We improve the methods employed for the calculation of Zp^ in the quenched 
theory [11] by including the mass term in the integrated Ward identity ()8.2j) . Of 
course, one still has to extrapolate to the chiral limit, but we will show that in 
practice the mass dependence of this normalization condition is extremely small. 

8.2.1 The vector current 

Our starting point is the continuum vector Ward identity (jS.lll and we take the 
region 7^ to be the entire Euclidean space for times smaller than some positive 
Xq. As internal operator we use the pseudo-scalar density at time zero and as 
external operator the pseudo-scalar density at time T > xq. The surface integral 
then becomes a spatial integral over the time component of the vector current. 
Applying (fTTl|l to (fT^ gives 

b\P\x) = -ie'"^P^{x) . (8.8) 

In this setup the vector Ward identity becomes 

((5tP^(0,u))P'^(T,v)) 
ie^'='^{^P'^(0,u)P"(T,v)^ 

The physical interpretation of this relation is that the isospin charge Jd'^xV^^(x) 
generates an infinitesimal isospin rotation of the state created by P'^, which trans- 
forms according to the vector representation of the exact isospin symmetry. 

If the open isospin indices on the right-hand side are contracted in a totally 
antisymmetric way with ie"''"^, the result is 

-2(^P"(T,v)P"(0,u)^ = ze'^'"=yd3x^P"(T,v)yJ'(xo,x)P"(0,u)^ . (8.9) 

On the lattice we now construct these matrix elements in the framework of the 
Schrodinger functional. We use the boundary field products ()5.13jl and ()5.14j) to 



- jd'^(v,\xo,^)P%0,u)P''{T,Y) 
/"d3x(p"(T,v)\/o'(xo,x)P^(0,u) 
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create initial and final states that transform according to the vector representa- 
tion of the exact isospin symmetry and insert the renormalized improved vector 
current ()4.H5j) . We demand that the correlation function 



/v (^o) = ^ E ^^^"'{0"'{VMx)0'^) , (8.10) 

X 

is equal to 

h = -^,{0'^0^) , (8.11) 

which corresponds to ()8.9|) when u and v are summed over all of space. In 
the improved theory it defines the renormalized vector current up to an O(a^) 
uncertainty 

f^ixo) = fi + Oia') , (8.12) 

since both correlation functions approach their (common) continuum limit with 
this rate. Note that we do not need to include the renormalization factors for the 
(multiplicatively renormalizable) boundary quark fields here because they appear 
on both sides and thus cancel. With the obvious definition of /v(a^o) 

/v(^o) = ^E^^'^''<^"(^i)o(^)^^) ' (8-13) 

X 

equation ()8.12jl then gives 

Zv(l + byam^)fl,{xo) ^ /i + 0{a') . (8.14) 

One can easily evaluate the contribution of the 0(a) counterterm appearing in 
the definition ()4.26|1 of the improved vector current to the correlation function 
/v(^o)- If we introduce the correlation function for the bare unimproved current, 



X 



it follows that 

Km 



oc 2^e°^"^C"^9^[^(x)(To^rV(a;)]C') . (8.16) 

X 

Since ao^ = f [7o,7^] is antisymmetric in its Dirac indices, we have (Too = 0. The 
contribution of the improvement term to the correlation function will therefore 
contain the expression 

5^4[^(x)ao,rV(a:)], (8.17) 
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which vanishes identically if we impose periodic boundary conditions in the spatial 
directions. We can therefore conclude that /v(a;o) = fvixo) and from equation 
it follows that 

ZyCgDil + 6vamq)/v(xo) = /i + 0{a^) . (8.18) 

By evaluating the correlation functions /i and /v(xo) through numerical simu- 
lation one is thus able to compute the normalization factor Zy{1 + hyarric^. In 
particular, to calculate Zy it suffices to consider the theory at vanishing quark 
mass. Equation ()8.18j) also implies that fyi^xo) is independent of xq up to cut- 
off effects. In a numerical evaluation one can therefore increase the statistical 
accuracy by averaging over a range of Xq values. More details will be given in 
Section 1831 

Derivation in the operator formalism 

The vector current normalization condition can also be derived in the operator 
formalism of the Schrodinger functional with boundary states S and S'. The 
fact that the charge is conserved implies that the corresponding charge operator, 
denoted by Q^, commutes with the transfer matrix T. The correlation function 
fl8.1()|) now reads 

= -^ie''''^{S'\d"'T^Q^d'\S) . 

Since Q'^lS) = we can replace Q^O'^ by its commutator and since the boundary 
fields are isospin vectors we have [Q^, O'^] = —i^bcdQd ^^^^ continuum. The phys- 
ical interpretation is that the isovector charge generates an infinitesimal isospin 
rotation as in ()8.8p . In an 0(a) improved lattice theory this argumentation holds 
up to O(a^) and therefore 

f^(xo) = —^e°^^{S'\d"'T^d''\S) + 0{a^) 

=2(5'"' 

= -^(OV) + 0(a2), (8.19) 

where in the last line we have switched back to the path integral form to recover 
the result (^n^ . 

8.2.2 The axial current 

We will start with a derivation of the normalization condition used in [11], where 
the quark is set to zero in the beginning. This will later be generalized to m 7^ 0. 
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Massless case 

For vanishing quark mass our starting point is the Ward identity ()8.3|) . Contract- 
ing the isospin indices gives 

= 2z(y;(y)axt) . (8.20) 

This relation is now transcribed to the 0(a) improved lattice theory with TZ 
being the region between the hyper-planes at Xo = ?/o±t and v = Q. With periodic 
boundary conditions in the spatial directions the surface integration in (|8.20p 
results in the difference between the axial charge at times ?/o =t t, i.e. 



\[(^R)o(^o+t,x) - (ylR)g(yo-t,x)J(AR)j5(y)axt 

X 

= 2i{{Vn)l{y)0,^,) + 0{a^) . (8.21) 



On-shell improvement is effective in ()8.21|) since the fields in the correlation func- 
tions are localized at non-zero distances from each other. Equation 1)8.211) is 
summed over y to obtain the axial charge. In the form of ()8.21|1 the integrated 
Ward identity requires three different time-slices, where the axial charge is in- 
serted. 

Using the conservation of the (renormalized) axial current the two insertions 
at yQ — t and ?/o, associated with the lower surface of 7^, can (simultaneously) be 
shifted to yo and yo+t. We thus arrive at 

a«^e'^''^((AR)«(x)(AR)^(l/)axt) = a'$^^((Vii)S(2/)axt) + 0(a2) , (8.22) 
x,y y 

where Xq = yo+t. Since t was arbitrary, the above equation holds for all insertion 
points xq and yo such that < yo < xq < T. 

At this point we have a relation that would allow us to calculate the axial 
current renormalization from the vector current renormalization. We can make 
things even simpler by choosing the field product O^xt such that the correlation 
function fyiyo) appears on the right-hand side of equation ()8.22|) . This amounts 
to setting 

axt = -^e^'^'C'^O'' , (8.23) 
since with this definition equation ()8.22|1 becomes 
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■^E^"''^'''(^"(^i^)o(^)(^r)o(?/)(^'^) (8-24) 

i3 



a 



6L6 ^ \ 
y 

,3 



/^(2/o) + 0(a2) . (8.25) 



The normalization condition for the vector current ()8.12|) allows us to replace 
fwivo) by /i. If we define the unrenormalized version of the correlator ()8.24|1 as 



x,y 



one can conclude from ()8.25p that 

Zlfij,{yo+t,yo)=h + 0{a') (8.27) 

for all times t > such that < yo and yo + t < T. As in ()8.12|1 the normalization 
of the boundary quark fields cancel. The axial current normalization constant 
Za can thus be determined by computing the two correlation functions /i and 
/aa(i/o+^5 l/o) at vanishing quark mass. 

Numerical simulations (see Section 18.3.111 show that the above relation has 
a very pronounced mass dependence. Since the quark mass was neglected from 
the very beginning in its derivation, this should not come as a surprise. In the 
following we will derive a normalization condition from the full PCAC relation. 



Non- vanishing PCAC mass 

Since we use a mass-independent renormalization scheme the normalization con- 
dition for the axial current can not be set up at finite quark mass. Instead our 
goal is to derive a normalization condition, which has a smaller mass dependence 
than ()8.27|1 . such that in practice the chiral extrapolation is easier. To this end we 
now use the axial Ward identity ()8.2j] with the same operator Oint = ^^(|/), where 
again the point y is somewhere in the interior of TZ. After anti-symmetrically 
contracting the isospin indices the result for an arbitrary mass m reads 

/ da^ix) e'''"^(A-{x)Al{y)0,^,) 
Jan ^ ' 

-2m^d^xe'^^^(p'^(x)4(?/)axt) = 2i(v,%y)0,^,) . (8.28) 
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Already here we note that the volume integral contains a contact-term, i.e. a 
contribution from the composite operator e"''^'^P"'{x)Al{y) as x approaches the 
(fixed) y. The operator product expansion implies that the leading contribution 
in the x^y limit comes from V^{x). From power counting it then follows that 
the composite operator cannot diverge faster than \x—y\~^, such that we receive 
a finite contribution under the four-dimensional integral over TZ. 

We choose the same region TZ and set i/ = as in the massless case. Together 
with an additional spatial integration over y the result is 

Jd'y J d^xe"^^ [(AS(i/o+t,x)A|;(|/o,y)axt) - (^^(yo-t, x)A|;(yo, y)axt)" 
-2m Jd^y y"d=^x^'"^*dxoe'^''^(p'^(xo,x)A|;(yo,y)axt) 

= 22 jd'y (y,%yo,y)0,^t) ■ (8.29) 

In the massless case the two contributions from the surface integration were the 
same due to current conservation and the anti-symmetric isospin structure. Here 
we have to use the partial conservation of the axial current to transform the 
second term of the surface integral. This will cancel the lower part {yo — t to yo) 
of the volume integral. The detailed calculation is given in Appendix O and the 
result is 

t Jd'y (Vo'(yo,y)axt) = 

Jd'y y"d=^xe"''^(Ag(2/o+t,x)4(yo,y)axt) 

-2m Jd'y y"d3x^'°^dxoe'^''^(p"(xo,x)A|;(2/o,y)axt) • (8.30) 

As before a normalization condition for the axial current on the lattice is obtained 
by demanding that eq. ()8.30j) in terms of the renormalized currents holds at 
non-zero lattice spacing. Inserting the same external operator ()8.23p as before 
will again allow us to replace the matrix-element of the vector current with the 
correlation function /i. As a result of the contact term 0(a) improvement fails 
in the correlator multiplying the mass term. On the lattice we therefore expect 
corrections of order am at finite mass in addition to the overall O(a^) uncertainty 
and thus have 



a 



x,y 

7 yo+t 



6L6 

xo=yo x,y 



/i + 0(a") + 0(am) . (8.31) 
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The weight factor w{xq) is needed to implement the trapezoidal rule in the dis- 
cretization of the time integral in ()8.H()j) . i.e. that the boundary terms Xo = yo 
and yo + t should only contribute with the weight 1/2. It is given by 

/ N J 1/2, Xo e {yo,yo+t} , . 

w{xo) = <( . (8.32) 

1, otherwise . 

If we set m = in ()8.3ip we immediately recover the normalization condition 
used in [11], i.e. ()8.25|) . With the definitions ()8.26|) and ()4.34|) the normalization 
condition can be written as 



^ 2mRa 



Zl{l + bAamqyfljs^{yo+t, yo) 

E ^(^o)$^e»''V'^^(o"^P^(x)(AR)^(l/)0^) = A + Oia') + 0{am) . 

xo=yo x,y 

Since there is no improvement term for the pseudo-scalar density we now define 
the new correlation function 

7 yo+t 

fkiyo+t,yo) = E w{xo)Y,e'^'^e^'^(o"P'^{x){A,)l{y)0^) . (8.33) 

xo=yo x,y 



From the PCAC relation we conclude that the product mP renormalizes with the 
same factor as A^. Thus, the final form for the normalization condition derived 
from the Ward identity with mass term is given by 

Zl{l + bAam^y{f{j^{yo+t,yo) - 2m/pA(2/o+t, yo) 



/i + 0(am) + 0(a2) . (8.34) 



8.3 Numerical computation 

Before going into the details of the simulations, we need to specify our choice of 
kinematic parameters in the numerical evaluation of the conditions ()8.18|) and 

The spatial boundary conditions are strictly periodic [6k = 0) and to ac- 
commodate the two insertion points in ()8.34|1 in a symmetric way, we choose a 
T = 9/4L geometry. The background field [Ck and C^) is set to zero and the 
improvement coefficients have the values specified in Sections 14.21 and 15.11 The 
axial current improvement coefficient ca is given by ()7.26|) . 

As first discussed in [11] and further detailed in [102], we need to evaluate the 
normalization conditions on a line of constant physics, keeping all length scales 
fixed. This ensures that the O(a^) ambiguities in the normalization factors vanish 
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smoothly when the perturbative regime is approached. In addition, the normal- 
ization conditions have to be set up at zero quark mass since we are implementing 
a mass-independent renormalization scheme as detailed in Section 14.31 

To keep the volume constant, we again employ the approximate scale setting 
given by cf. eq. ()7.2|) . We point out that in the present deviation 
from the line of constant physics can change the result by O(a^) only, whereas 
in the computation of ca such deviations would show up at 0(a). In addition to 
three lattice resolutions matched in this way, we simulated at three larger values 
of (5 and fixed L/a = 8, which results in very small volumes. This was done 
in order to verify that our non-perturbative estimate smoothly connects to the 
perturbative predictions, eqs. ()4.36|) and ()4.37|) . 

Our "reference" volume is an 8'^ x 18 lattice at /5 = 5.2, which corresponds to 
the coarsest lattice spacing. At all other lattice spacings the systematic error in 
Za and Zy due to a mismatch in the volume is estimated by varying the lattice 
resolution L/a. To check for smoothness of Za((7o) additional simulations were 
performed in an unmatched 8'^ x 18 volume at (3 = 5.29 without an estimate of 
the systematic error. 

These last simulations as well as all but the heaviest run a.t (3 = 5.2 were 
done with the PHMC algorithm. All others employ the HMC with two pseudo- 
fermions as discussed in Section 16.11 For all runs it could be verified that the 
problems described in Sections 16.51 and are absent. 

For reference we collect again the precise non-perturbative definitions of Zy 
and Za- In the simulations both the PCAC mass and the estimate for Zy are 
averaged over a few time-slices in the middle of the lattice. No such time average 
is performed for the correlation functions /^^ and /p^. 



Za(^o) = lim X — , (8.35) 

V /1a(2T/3, r/3) - 2mfU2T/3, T/3) 

^v(^?o) = limi^X^y^, (8.36) 



m 



I Jl^ do\fA{xo)+gA{T-xo) +acAdod^ fp{xo)+gpiT-xo) 



Nt 

X0=tl 



do 


fA{xo)+gA{T-xo) 


+ ac Adod^ 


fp{xo)+gp{T-xo) 


2 


fpixo)+9piT-xo) 





^he.e ('■ = X''' = S'^' = ' (8.38) 

y ti = ^, t2 = ^, Nt = 4: for T odd. ^ ^ 

8.3.1 Implementation notes and quenched example 

The starting point in the implementation of the correlation functions for the axial 
normalization condition were the APEmille TAO codes of the ALPHA collabo- 
ration. These already include the correlation functions /a, /p, Qa and gp. 
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am 

Figure 8.1: Mass dependence of the old and new normalization conditions for 
the axial current evaluated in the quenched approximation. The parameters are 
/? = 8.0, volume = 8^x18 and n e [0.1325,0.133168]. 



The correlation functions /v, /aa ^ind /p^ have been implemented in the 
quenched code (which uses SSOR preconditioning [118]) as well as the even-odd 
preconditioned HMC and PHMC codes. An additional version (derived from the 
quenched one) exists, which measures the correlation functions on a saved gauge 
configuration. A detailed discussion of all correlation functions and in particular 
of ^'^d /pA and their Wick contractions is given in Appendix IB. II 

The code for /aa as well as the analysis program have been checked against 
those used in [11]. Agreement was found both statistically (i.e. at a given set of 
bare parameters) as well as on a given gauge configuration. 

As a first test of the new axial current normalization condition, we consider 
sets of quenched 8^x18 gauge configurations at /3 = 8.0 and various values of the 
hopping parameter n. The lattice spacing and hence the volume are extremely 
smalP such that the fluctuations of the gauge field are strongly suppressed and 
one can (in the Schrodinger functional) easily simulate at the critical point. 

The results are shown in Fig. 18.11 where one can see that the lightest mass 
(with K = 0.133168) is exactly at the critical point. We show the estimate for 
Za defined in fl8.35|) as well as the one from the old condition [11], where the 
term mfl,^ is neglected. As expected from the argumentation in the previous 
section, the new condition has a significantly smaller mass dependence. At these 

^In the quenched approximation /3 = 6.0 implies a ~ 0.1 fm and hence corresponds to (3 — 5.2 
in the two-flavor case. 
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parameter values using the new condition reduces the linear coefficient of a fit 
in am from —14.9(8) to 0.7(2), which is consistent with our 0(am) expectation. 
With a value of 0.97(2) the slope in the estimate of from (jS.Hfjj) . which is not 
shown in the plot, is similar to that of with the new condition. 

Due to the large slope, an extrapolation of results from the old normalization 
condition of the axial current should take into account the error in the quark mass 
as well as the correlation between am and the Za estimate. In contrast, with 
the small slope from the new condition the extra- or interpolation is essentially 
fiat and uncertainties in the location of the critical point do not propagate to the 
determination of Z^. 

Note that in the shown quark mass range the old normalization condition al- 
ready shows non-linear effects in am. Since a is extremely small here, the physical 
quark mass is already quite large and such a behavior is thus not unexpected. 

8.3.2 Results for the normalization factors 

The results of all simulations are collected in Appendix^ while in Table IHTTl we 
show only the results of the chiral extrapolations ()8.35|) and ()8.36p . 

In most cases at least one simulation is very close to the critical point and in 
some cases we actually perform an mterpolation in am. As can be seen from the 
table, this results in a very precise estimate of the critical hopping parameter Kc. 
The latter is obtained by extrapolating the PCAC mass am linearly in 1/k. 





L/a 


T/a 




Za 


Zy 


5.200 


8 


18 


0.135856(18) 


0.7141(123) 


0.7397(12) 


5.500 


12 


27 


0.136733(8) 


0.7882(35)(39) 


0.7645(22)(18) 


5.715 


16 


36 


0.136688(11) 


0.8037(38) (7) 


0.7801(15)(27) 


5.290 


8 


18 


0.136310(22) 


0.7532(79) 


0.7501(13) 


7.200 


8 


18 


0.134220(21) 


0.8702(16)(7) 


0.8563(5)(45) 


8.400 


8 


18 


0.132584(7) 


0.8991(25)(7) 


0.8838(13)(45) 


9.600 


8 


18 


0.131405(3) 


0.9132(11)(7) 


0.9038(3)(45) 



Table 8.1: Results of the chiral extrapolations to extract ()8.35|) and ()8.36|) . 



As already mentioned the volumes of the simulations at /3 = 5.2, 5.5 and 
5.715 are matched in the same way as for the ca simulations, cf. eq. ()7.2j) . while 
for (3 > 7.2 the volumes are much smaller.^ The second set of errors for Za/v 
represent an estimate of the systematic uncertainties due a mismatch in the lattice 

^The (compared to ChapterEj) different f3, such that t(/3) — 0.5 is due to the use of a version 
of eq. (j7.2(l . which differs by 0(5q) terms. 
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spacing/volume. They will be discussed in Section [8.3.31 An exception are the 
simulations at /5 = 5.2, which is the reference volume and (3 = 5.29, which was 
performed only to qualitatively confirm the observed rapid change of Zx^Qq) in 
this region of the coupling. 

Fig. 18.21 shows the critical mass arrtc = (2/tc)~^ — 4 as determined in the sim- 
ulations and compares it to the (continuum) perturbative prediction. The non- 
monotonic behavior of the additive mass renormalization at Qq > 1 is similar to 
the one found in the quenched case (see e.g. [11]) and is clearly a non-perturbative 
phenomenon. In contrast, for the largest values of (3 simulated, the numerical 
results are already close to the the two-loop formula ()4.12|) . 

The chiral extrapolation of our estimate for Za at the coarsest lattice spacing, 
(3 = 5.2, is plotted in Fig. 18.31 Apart from the significantly larger errors the 
situation is very similar to the quenched one from Fig. 18.11 showing that also 
in the dynamical case the new normalization condition has only a very small 
mass dependence. As a result, the errors on am and the correlation of am with 
the estimate for Za are irrelevant for the extrapolation. Note however, that the 
mass range in lattice units is the same in both plots, which implies that (ignoring 
renormalization factors of order one) the physical masses in the dynamical case 
are much smaller since a is much larger. As a result, also for the old normalization 
condition all mass effects now show a linear behavior. 

The slopes of the two estimate as functions of am are also very similar to the 
quenched case. Due to the large statistical errors the slope for the new condition 
is consistent with zero (0.4 with an error of 1.2) and for the old condition we 
obtain -14(1). Again, the slope of the estimate for Zy (0.64(9), not shown) is 
comparable to that from the new condition for Z^- In the quenched case [11] the 
slope (in am^) of fi/{Zyfy) could be used to determine the O(amq) improvement 
coefficient by, cf. eq. ()8.18|1 . This was possible because for A^f = the bare 
and modified gauge couplings coincide {bg = 0). This is no longer true in the 
dynamical case [63,71] and hence the slope one finds is not entirely due to by. 

In order to keep the discussion transparent, we will now restrict ourselves to 
the new normalization condition. The final results for Z^ and Zy are shown in 
Fig. 18. 41 as a function of Qq. The errors plotted there are obtained by linearly sum- 
ming the statistical and systematic errors. Also shown are the 1-loop estimates 
flOF)|l and (jOTjl . 

One can see that our results for both Z^ and Zy lie on smooth curves and 
we can therefore (as for ca) describe the non-perturbative data in form of a 
rational function of Qq (Fade fit). In this fit the difference of the linear terms of 
numerator and denominator is constrained to the 1-loop value, which ensures a 
smooth matching with perturbation theory as Qq approaches zero. 

The fit is performed with the sum of the errors from Table 18.11 and the data 
points aX (3 = 5.29 are excluded since no estimate of the systematics is available 
in this case. With the perturbative input ()4.36|1 and ()4.37|1 we thus arrive at the 
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Figure 8.2: The critical mass arric (errors are much smaller than the symbols) as 
a function of the bare gauge coupling Qq. For < 0.8 the non-perturbative data 
start to approach the two-loop prediction (|4.12|) from [55]. 
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Figure 8.3: The chiral extrapolation of the result from the axial current normal- 
ization conditions at (3 = 5.2. 
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Figure 8.4: The normalization factors Zx and Zy as a function of g^, together 
with the 1-loop results (dashed line) and the Pade interpolating formulae (solid 
line). The data points at /3 = 5.29 (lighter color) are excluded from the fit. 



interpolating formulae 

= l-0.918«g + 0.062,j + 0.020 ,g 

1-0.8015^2 ' V ; 

7 (r?\ l-0.6715^g + 0.0388^0^ 

= 1 - 0.5421^0^ ■ (8-40) 

The low value of Za at /3 = 5.2 requires that the Za data are fitted with a 
polynomial of third degree in the numerator, while the Zy data can be well 
described including only quadratic terms in in the numerator. For Zy we 
find agreement at the 1% level with the results from [119], where isospin charge 
conservation is imposed for nucleon matrix elements of the local vector current 
in large volume. 

The maximal absolute deviations of the fit from the simulation data are 0.0056 
(corresponding to 0.7 a or 0.7%) for Za at /5 = 5.5 and 0.004 (0.9 a or 0.5%) for Zy 
aX (3 = 5.715. Considering that all the data in one fit are statistically independent 
since they come from different simulations, these fluctuations are rather small. 
For future application we propose to either use the data from Table 18.11 directly 
or the interpolating formulae ()8.H9j) and ()8.40j) . To Zy we ascribe an absolute 
error of 0.005 and for Za the absolute error decreases from 0.01 at /5 = 5.2 to 
0.005 at /3 = 5.7. 

We note that for a lattice spacing of roughly 0.1 fm, corresponding to [3 = 5.2, 
our non-perturbative estimate of Za is almost 20% smaller than the 1-loop value. 
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while in the quenched case [11] this difference was only 10% at the same lattice 
spacing.^ 



8.3.3 Systematic uncertainties 

Close to the continuum the dependence of the normalization factors on the lattice 
size is expected to be of order (a/L)^ [11] in the improved theory. This implies 
that effects in Za and due to deviations from the line of constant physics, i.e. 
the constant physical volume condition, should be strongly suppressed. 

To check for these effects, a.t (3 = 5.5 and 5.715 the simulations closest to the 
critical point were repeated on smaller lattices {L/a = 8 instead of 12 at /3 = 5.5 
and L/a = 12 instead of 16 at /3 = 5.715). From our approximate matching (see 
Chapter [7j) we estimate the uncertainty in L (measured in units of L at (3 = 5.2) 
to increase up to at most 10% in the range 5.2 < (3 < 5.715. We assume that 
the uncertainty in L grows linearly in j3 and therefore assign a 6% error to L at 
[3 = 5.5 and 10% at /5 = 5.715. Together with the difference in the estimates of 
Za and Zy this gives the systematic errors quoted in Table 18.11 through a linear 
propagation of the error. 

Our simulations confirm the expectation of small volume dependence. The 
largest difference is seen at /3 = 5.5 {L/a = 12 vs. 8), where it is 0.022(7) 
corresponding to a 3% effect in Za for a 33% change in L. 

For the runs with (3 > 7.2 the matched L/a would be extremely large. On 
the other hand the volume dependence of the normalization factors should be- 
come smaller as we are approaching the perturbative regime. In practice the 
simulations are thus performed at L/a = 8 and the systematic error is estimated 
from additional runs at the coarsest of these lattice spacings, i.e. at /5 = 7.2, 
by taking the difference of Z\ and Zy between L/a = 8 and 16. While for the 
Za estimate the volume dependence is hardly visible, in the case of Zy the large 
volume (L/a = 16) results in a (statistically) significantly lower value. This is 
the reason for the larger systematic error of Zy for j3 > 7.2. 

As already mentioned, the simulations aX (3 = 5.29 (again with a 8'^x 18 lattice) 
were done to check that the large difference in Z\ when going from j3 = 5.2 to 5.5 
is bridged smoothly. At this value of (3 the matched volume would be L/a = 9, 
which is not accessible with our machine (APE geometry restrictions) and code 
(even-odd preconditioning). Although the results were excluded from the fits 
for this reason, also the (3 = 5.29 data are well reproduced by the interpolating 
formulae (see Fig. I8.4|l . which implies that possible systematic effects are not 
visible at the given statistical precision. 



"'Although m [11] Z™" (see next section) was considered, this comparison is still meaningful 
since the difference from our definition of was found to be negligible already for a ~ 0.1 fm. 
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Figure 8.5: Quark diagrams for with non-vanishing isospin factors. The gray 
diamonds indicate the insertions at times yo and xq. 



8.3.4 Comparison with alternative normalization 
conditions 

In this section we consider two variants of the definition of fl8.35|) and study 
the difference in the results for Za(5'o)- This provides us with an estimate for the 
magnitude of the cutoff effects to be expected at the considered lattice spacings. 

The Wick contractions of /aaI^Oj^o) ^-re derived in Appendix IB. II with the 
result that for the chosen isospin assignment only the quark diagrams shown 
in Fig. 18.51 contribute. Diagrams related by an exchange yo Xo have the 
same isospin factors with opposite signs. This follows directly from ()8.2(ij) . where 
such an exchange corresponds to e"'^'^*-* e^"'^. Among the quark diagrams are the 
two disconnected diagrams b) and c), where no quark propagator connects the 
boundaries. 

Following [120], we now argue that in the massless continuum limit the con- 
tribution of the disconnected diagrams vanishes. Thus we consider diagram b) 
without improvement, i.e. the current insertions are just axial charges. 

In this diagram, termed [gf]AA{xo,yo) in Appendix IB. 21 the spatial insertion 
points X and y are integrated over such that it is a function of Xq and yo only. 
In fact, since we are in the chiral limit, the axial charge is conserved and hence 
the diagram is independent of the insertion points in the two regions xq < yo and 
^0 > Vo- If fhe two insertion meet, contact terms may arise, which we need to 
treat separately. 

To this end we restrict the spatial integration to |x — y| > e and let the 
insertion times approach one another from either region. No contact terms can 
appear due to the finite spatial separation. If the integrand has a divergence 
weaker than |x — y|~^, the contribution to the spatial integral from the region 
|x — y| < e vanishes as we make it smaller and smaller. In this case we can 
take the limit e ^ and conclude that the order of xo and yo does not play 
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any role. This would imply [gf]AA{xo,yo) = [gf]AA{yQ,xo), i.e. the diagrams b) 
and c) have the same value. Since their isospin factors have opposite signs, their 
contribution to /aa cancels. 

This argumentation shows that in the massless continuum theory the dis- 
connected quark diagrams b) and c) do not contribute if the composite field 
74o(a;)y4o(y) has a divergence weaker than |x — y|~^. To decide this we consider 
the matrix elements of two axial charges between pseudo-scalar states and assign 
the flavor quantum numbers such that diagram b) is the only Wick contraction. 
Using four valence quark species u, d, s and c this amounts to 



It is now immediately clear that such a flavor assignment excludes a single 
quark bilinear as the leading contribution in the operator product expansion 
of Ao{x)Ao{y). Hence the latter has (if any) a divergence weaker than |x — y|~^ 
in the limit |x — y| e and the contribution from the excluded integration region 
vanishes. 

Since the correlation functions approach their continuum value with a rate 
proportional to a^, we can conclude that on the lattice the contribution from the 
disconnected diagrams is a cutoff effect of this order. 

If, in contrast, we would consider the diagrams a) and d), which also differ 
by a sign, they are still independent of the insertion times in the regions xq < yo 
and Xq > yt). However, in the above argumentation the flavor structure would be 



and hence the operator product could mix with e.g. (Vo)dc or the scalar density. 
This, by dimensional analysis, produces a factor |x — y|~^ in the integrand at 
leading order of the OPE and thus, after integrating, a logarithmically divergent 
contact term. Therefore, no statement about the contribution of the connected 
diagrams can be made from this argument. 

In more physical terms this can be understood directly from the quark dia- 
grams. Even if the insertion points coincide when we move them around (on a 
fixed gauge background), in b) there is no propagator connecting them and the 
axial charges thus don't "see" each other. Doing this with diagram d) changes 
the propagator J^x y ^i^'i v) ^i^d in particular it gives a contact term at zero 
separation. 

As is shown in Fig. 18.61 we can numerically confirm that the contribution of 
the disconnected diagrams to vanishes in the continuum limit. According to 
Appendix IB. 21 the contribution is 







disconnected 



[fi'/]AA(xo, yo) - [5'/]aa(2/o, a;o) 



+acA<9o i - [fl'/]pA(a;o,yo) - [fi'/]Ap(2/o, a;o) 
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Figure 8.6: Data from the simulations at /3 = 5.2, 5.5 and 5.715 extrapolated 
to the chiral limit. As expected the contributions from the disconnected diagrams 
vanishes in the continuum limit with a rate that is faster than linear in a. 



+acyi(9^ |[5f/]Ap(xo,yo) + [gf]pA{yo,xo)^ 
WcldS ~dl { [(7/]pp(yo, xo) - [^7/]pp(xo, 2/o) \ , (8.43) 



where we have Uq = T/3 and Xq = 2T/3. In Fig. 18.61 we plot data from the 
matched simulations and neglect any systematic effects from volume mismatch. 
One can clearly see that the above term approaches zero faster than linear in the 
lattice spacing and is in fact compatible with an effect. 

The correlation function multiplying the mass term, /p^, only influences the 
slope of the chiral extrapolation but not the result for Za in the chiral limit. 
We are therefore free to also drop the disconnected diagrams for this correlator. 
An estimate of Z\ using ()8.35j) . where only the connected part of the correlation 
functions enter, should therefore agree with the original definition up to O(a^). 
We will denote this by 

In the quenched case [11] was also determined from the connected quark 
diagrams only. The authors approach the above conclusion from a different point. 
Their argument exploits the freedom to choose different external operators in the 
Ward identity (jH221)- Instead of ^2^, one could use 

(^ext = ^,0[liT%0, + ^,JC[l{T%}C, , (8.44) 
where 0', = a'J2 Ci^H^'i^) ' 0, = a'J2 ^HlM^) , (8-45) 

u,v u,v 

= a^J^f (u)75C'(v) , /C. = a^^C(u)75e(v) • (8.46) 



u,v 
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Figure 8.7: Changing the definition of Za at O(a^) has significant effect on the 
result for bare gauge couphngs ^ l-l- 



Here i and j are flavor indices and ^ denotes the boundary field for a third quark 
species, which is taken to be an isospin singlet. Therefore no Wick contraction 
exists that connects it to the axial current with the result that no disconnected 
diagrams can appear. Since the external operator in the Ward identity is arbi- 
trary, the resulting Za agrees with the original one up to cutoff effects and one 
thus again reaches the same conclusion. Note that the first argument takes place 
at the level of quark diagrams and is thus valid for any number of valence and 
sea quarks, whereas the external operator ()8.44j) explicitly assumes the existence 
of three valence quarks flavors. 

If the improved axial current is inserted into the correlator /^^, we obtain 
expressions proportional to c\ as in the last line of ()8.43p . Since we implement 
improvement only to 0(a), we could have decided to drop these as just another 
source of O(a^) ambiguities. Thus, keeping only improvement terms linear in ca 
provides us with another estimate of the axial current normalization constant, 
which we denote by Z^™.^ 

We can now repeat the whole analysis - chiral limit at each value of (3 and 
estimate of volume dependence - to extract Za^Qq) from the estimates Z'^^ and 
Z^^. The result of this exercise is shown in Fig. 18 .71 where it is compared to the 
interpolating formula ()8.H9j) for the original estimate. We note that for (3 > 5.7 
[qq < 1.1) the central values agree with the interpolating formula ()8.H9|) . The 

^One might be tempted to also consider but in practice the contribution of the c\ 

term to the connected part is negligible and hence ~ 
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increased errors for at /5 > 7.2 are due to tlie ratlier large dependence on 
L/a we found at /? = 7.2 (see Appendix lD|l . 

However, much more noticeable are the O(a^) effects in the range oi Qq ~ 1.1 
and above, corresponding to /3 < 5.4. The agreement with the original determi- 
nation in form of the Pade fit is rapidly lost and ai (3 = 5.2 we see a (statistically 
significant) deviation of 15% (for Z'^'^). 

This might be surprising, given that in Fig. 18.61 the /? = 5.2 estimate of the 
disconnected part is only two standard deviations from zero. It is, however, of 
similar magnitude as the connected part and the small error in the resulting Za 
is due to the correlation of numerator and denominator in ()8.35|1 . 

From the above observations and Fig. 18.71 it is clear that it is essentially the 
last term in ()8.4Hj) . which is responsible for the monotonic behavior in our original 
estimate. The large contribution of this term is in turn related to the fact that in 
the volumes and time extensions we consider here, the /p correlator still shows a 
strong exponential decay even at the critical point. This means that it receives 
large contributions from excited states in the pseudo-scalar channel [114] and 
therefore its effective mass large. 

These deviations and even the apparent non-monotonic behavior one obtains 
from Z^'^ and Z^^ do not signify any theoretical problems. It is, however, another 
sign that even after a successful implementation of the improvement programme 
the lattice artifacts in the two-flavor theory can be large for /3 < 5.4 and a 
continuum extrapolation remains necessary to make reliable physical predictions. 
In particular, care should be taken that such extrapolations are not dominated 
by the data from simulations with f3 = 5.3 or lower. 

In the quenched case the O(a^) effects from the disconnected diagrams are 
significantly smaller at a lattice spacing a ^ 0.1 fm. In fact, at a statistical 
precision of the order of 1% these effects are not visible at all. 

8.4 Summary 

We have formulated an integrated Ward identity, derived from the isovector chiral 
symmetry of the continuum action, as a normalization condition for the axial 
current on the lattice. Through the use of the massive Ward identity we were 
able to improve the method from [11] with respect to its mass dependence. 

The normalization condition was implemented in terms of correlation func- 
tions in the Schrodinger functional framework and evaluated on a line of constant 
physics in order to achieve a smooth disappearance of the O(a^) uncertainties. 
Through additional simulations at very small lattice spacings and volumes we 
verified that our non-perturbative definition approaches the perturbative predic- 
tion at small bare gauge coupling. Simulations were done at or near the critical 
point and owing to the new normalization condition, which keeps track of the 
mass term in the PCAC relation, any chiral extrapolations were extremely fiat. 
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Systematic effects due to deviations from the constant volume condition are also 
estimated and turn out to be small. The results are well described by an inter- 
polating formula ^A(fi'o) range of bare couplings considered. 

Enforcing isospin symmetry with the same programme, we obtain at the same 
time a non-perturbative determination of the normalization constant Zylg^) of 
the local vector current. 

We found large O(a^) uncertainties in at /5 < 5.4 by varying the definition 
of Za, which indicate that despite improvement the cutoff effects in this range 
of lattice spacings can be large. Together with the algorithmic issues discussed 
in Chapter IHl these findings corroborate the worries expressed in [85] about the 
status of simulations with improved dynamical Wilson fermions at the currently 
accessible lattice spacings. This merely emphasizes that cutoff effects in physical 
quantities can be controlled only if a continuum extrapolation with several lattice 
resolutions is performed. 

Finally, the method employed here can also be used to obtain Za with other 
actions formulated in the Schrodinger functional. This includes in particular 
the 0(a) improved three flavor theory with either Iwasaki or plaquette gauge 
action [60,116,117]. 



Chapter 9 
Conclusions 



The main result of this work is the non-perturbative determination of the axial 
current improvement and normalization constants for the lattice theory with two 
flavors of dynamical Wilson fermions. 

With this regularization we have thus shown that normalization conditions can 
be imposed at the non-perturbative level such that isovector chiral symmetries 
are realized in the continuum limit. Since we work with an improved theory, 
chiral Ward-Takahashi identities are then satisfied up to O(a^) at finite lattice 
spacing. 

In the course of the implementation of this programme algorithmic difficulties 
were encountered in the numerical simulations. A more detailed investigation 
revealed that those are related to features of the infrared spectrum of the Wilson- 
Dirac operator. More precisely, it was found that at a given lattice spacing the 
inclusion of the fermion determinant enhances the probability of finding very 
small eigenvalues, which is contrary to the usual intuition. This effect disappears 
rapidly as the lattice spacing is decreased and is thus interpreted as a lattice 
artifact. 

This spectral distortion is responsible for problems in the molecular dynamics 
evolution of the hybrid Monte Carlo algorithm as well as in the efficient sampling 
of observables sensitive to small eigenvalues. One of the important results is that 
in this situation it can be advantageous to deviate from importance sampling. In 
practice we found the polynomial hybrid Monte Carlo algorithm with reweighting 
to be an efficient tool in the parameter range we considered. 

The non-perturbative determination of Cs^, the improvement coefficient of 
the Wilson fermion action, has been known for some time already [59]. However, 
for a complete removal of the 0(a) lattice artifacts also the composite fields 
need to be improved. We have completed this programme for the axial current, 
accomplishing a necessary step on the way to a non-perturbative normalization 
of the latter. 

The non-perturbative determination of the axial current improvement con- 
stant ca was implemented in the Schrodinger functional framework using spatial 
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wave functions in the composite boundary fields. In this way we could obtain a 
large sensitivity by explicitly projecting to the ground and first excited state in 
the pseudo-scalar channel and requiring the current quark mass derived from the 
PCAC relation to be the same in both cases. With this method the contribution 
of excited states is well controlled and we were thus able to ensure that the un- 
avoidable 0(a) ambiguities in ca remain small. To achieve a smooth disappearing 
of the latter in the perturbative limit we evaluated the improvement condition on 
a line of constant physics. Lacking a non-perturbative estimate of the evolution 
of the lattice spacing with the bare gauge coupling, we used (in a small range of 
bare gauge couplings) an estimate based on the three-loop lattice /3-function for 
improved Wilson fermions. 

In the next step we implemented the non-perturbative renormalization of 
the improved axial current. Similarly to the determination of the improvement 
coefficient this is also achieved by enforcing a continuum Ward identity at finite 
lattice spacing. We improved the method used in the quenched case [11] to obtain 
a normalization condition whose dependence on the quark mass is only a cutoff 
effect, thus alleviating the need for a chiral extrapolation. 

Similar to ca, beyond perturbation theory Z^igo) is affected by an intrinsic 
O(a^) ambiguity, which vanishes smoothly in the continuum limit if the normal- 
ization condition is evaluated on a line of constant physics. However, in the case 
of Z\ we need to evaluate this condition at zero quark mass at each lattice spac- 
ing and Za is thus the result of an extra- or interpolation in the quark mass. To 
match the volumes we use the same estimate of da/dg^ as in the determination 
of Ca and with additional simulations at large (3 we verify that our definition of 
Za smoothly matches with perturbation theory. With little additional effort we 
also extracted the local vector current normalization factor Zylg^) from the same 
simulation data by enforcing isospin symmetry at finite lattice spacing. 

The dependence of a non-perturbative definition of Za on the lattice size is 
expected to be of order (a/L)^ close to the continuum. In our estimate of the 
systematic effects in Za due to deviations from the constant physics condition, we 
indeed observed only very little volume dependence. However, at lattice spacings 
a > 0.07 fm large ambiguities in Za were found by changing the normalization 
condition at O(a^). Together with the cutoff effects in the spectrum (Chapter IH]) 
these findings add to the evidence collected in [85] that the lattice artifacts with 
dynamical improved Wilson fermions at the current lattice spacings are much 
larger than initially expected from quenched experience. One might even sus- 
pect the proximity of a phase transition in bare parameter space in analogy 
to [84]. Such an extreme example of a lattice artifact would surely invalidate the 
Symanzik expansion ()H.17|1 . 

Setting these worries aside, the results for ca and Za can now readily be ap- 
plied. The most immediate application would be to 0(a) improve and renormalize 
the bare pseudo-scalar decay constants computed in [97,98,115]. As already men- 
tioned in Section Ell using the non-perturbative ca lowers the result from [115] 
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by 10% at /3 = 5.2. With a 20% decrease compared to tlie 1-loop estimate (again 
at /5 = 5.2) the non-perturbative Za will also have a large effect when used to 
renormalize existing data. 

The methods developed and implemented here can easily be applied to other 
actions formulated in the Schrodinger functional framework. This includes im- 
proved gauge actions [28] as well as theories with more than two dynamical quark 
flavors [60,116,117]. 

Within the research programme of the ALPHA collaboration, the results ob- 
tained here are an essential step in the computation of the pseudo-scalar meson 
decay constant Fps needed to reliably convert the A parameter from [2] into phys- 
ical units. In the short term, together with data from [14] Za(5'o) used in 
a fully non-perturbative calculation of the strange quark mass [15] following the 
strategy of [13, 72]. 
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Appendix A 

Improved action for the 
Schrodinger functional 



The 0(a) improved action of the QCD Schrodinger functional [63] is given by 
(the functional dependence on U, ip and ip is suppressed) 

Simpr = Sg + Sf + 6S,a + 6Sf^b ■ (A.l) 

The improvement of the gauge action ()3.4|1 is implemented through 

S,[U] = \Y.'^ip)TT{l-U,} , (A.2) 

where the weight for the oriented plaquette p is given by 

{1, for p in the interior of the SF 
Ct, for p touching the xq = or T boundaries (A.3) 
|cs, for p in the boundary. 

A sufficient condition for the Cg term not to contribute are Abelian and spatially 
homogeneous boundary fields Ck and C^. This condition is not only fulfilled for 
vanishing background field (the setup employed here), but for all boundary gauge 
fields so far considered in the SF [2,59]. 

The fermionic action is given by ()4.4|) in the interior of the SF cylinder and 
the boundary improvement terms ()5.9|1 and ()5.10j) are transcribed on the lattice 
as 

ssf,, = (£3-i)[a(x)+a;(x)] + (Q-i)[a,(x)-a^(x)], (A.4) 

X 

with the boundary and near-boundary operators 

a(x) = p(x)7fcVfcp(x) , (A.5) 



105 



106 



Improved action for the Schrodinger functional 



a:(x) = p'(x)7.V,p'(x) , (A.6) 
ai(x) = |^(x)P+VS^/'(x) + ^(x)V^P_V'(x)| , (A.7) 

(5;(x) = |^(x)P_Vo^(x)+^(x)VoP+^(x)) . (A.8) 

L J xo=T-a 

While Ot and 0[ couple to the dynamical variables close to the boundaries, Og 
and O'g depend on the fermionic boundary fields p, p, p' and p only. Since they 
are set to zero, a knowledge of Cg is not necessary. 

The entire fermionic part of ^A.l\ is expressed in {5D summarizes the volume 
and boundary improvement) 

5/,impr = Y,^^^)[D + 5D + mo] V'(^) , (A.9) 



X 



and for a given gauge field configuration ?7^(x) the quark propagator S{x] y) is 
thus defined as [61] 

{D + 5D + mo)S{x; y) = a'^S^y , < a;o < T , (A.IO) 

and the boundary conditions 

P+S{x-y) =P-S{x-y) =0, (A.ll) 



xo=0 



xo=T 



where color and spin indices have been suppressed. The basic two-point correla- 
tions are obtained [61] by differentiating the generating functional with respect 
to the source fields. The result can be expressed in terms of the quark propagator 
S{x;y). 



[tp{x)ij{y)]p 


= S{x]y) 


(A.12) 


[V^(^)C(y)]p 


= dtSix;a,y)UoiO,yr'P+ 


(A.13) 


[^(:^)C'(y)]p 


= CtS{x]T-a,y)Uo{T-a,y)P- 


(A.14) 


[C(x)^(2/)]p 


= CtP_f/o(0,x)S'(a,x; y) 


(A.15) 


[C'(x)^(2/)]p 


= CtP+Uo{T-a,x.y^S{T-a,:>L;y) 


(A.16) 


[C(x)C'(y)]p 


= £?P„f/o(0, x)5(a, x; T-a, y)f/o(T-a, y)P_ 


(A.17) 


[C'(x)C(y)]p 


= 5?p^f/o(T-a,x)-i5(T-a,x;a,y)f/o(0,y)-ip+ 


(A.18) 



The correlators ^ and [('('] p are not needed since the isospin factor for such 
a Wick contraction vanishes for isovector boundary field of the type ()5.1Hj) and 



Appendix B 
Correlation functions 



B.l Summed two-point correlators 

We introduce the quark propagator from the lower boundary to the point x 

^(^) = ^E^(^5«'y)^o(0,y)"^P+ (B.l^ 



? ^3 



y 



y 

^ Yl P+UoiO, y)75^(a, y; xh, , (B.2) 



y 

the quark propagator from the upper boundary to the point x 

^(^) = ^Y.^{x-^-(^.y)V^{'r-a,y)P. (B.3) 
y 



^ ^(^y = |0^$^i^-f/o(T-a,y)-S55(T-a,y;x)75, (B.4) 

y 

and the boundary-to-boundary propagator 

St = Cta2^-f/o(T-a,x)-ip+;S(r-a,x) (B.5) 

X 

^ = gta^^-5(T-a,x)^P+[/o(T-a,x) . (B.6) 

X 

Of course St can also be obtained through a summation of R{x) 
IbTI c^a^ 



755t75 = -^5y^^-f/o(T-a,x)-ip„755(T-a,x;a,y)75f/o(0,y)-ip_ 



x,y 

W c,a^5^-P(a,y)tf/o(0,y)-ip_^ 
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Correlation functions 



With these definitions we can express the basic two-point functions as 



a' 
L372 



y y 



(B.7) 



L3/2 



E [^(^)C'(y)]p'^ 0E^(^'^-«'y)f^o(T-a,y)P_,W:R(a:) . (B.8) 



The >S'(x) and R{x) propagators inherit the 75-Hermiticity from S{x] y) 



a 



x] 



a^ct 
L3/2 

a^ct 
L3/2 



J2[s{a,y;xyUo{0,yr'P- 

y 

E [75^(0, y;a;)t75f/o(0,y)-ip+75]^ 



a 



E[C'(y)^(^)], 



755(x)^75 , 

2_^P+Uo{T-a,y) S{T-a,y; 



(B.9) 



L3/2 

L3/2 



X) 



J2[s{T-a,y;xyUo{T-a,y)P^ 

y 

3 ~ 

J^Yl [^5S{T-a,y;xy-f5Uo{T-a,y)P. 



75 



75 



75i?(a;)S . (B.IO) 
For the spatially averaged boundary-to-boundary correlators one obtains 

^E [C(y)C'(x)] f$ E y' ^)f^o(T-a, x)P_ 



L3/2 



x,y 



x,y 



a6£2 



L3/2 



I E^s^+^o*^*^' y)^55'(a, y; T-a, x)75?7o(T-a, x)P+75 



x,y 



,3- 



«"Ct E 75'S'(T-a, x)^t/o(T-a, x)P+75 

X 



E6J 

= -75'3r75 



(B.ll) 



B.l Summed two-point correlators 
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%I1 [C'(x)C(y)] f$ P^U,{T-a, ^)-'S{T-a, x; a, y)f/o(0, y)-^P, 



13/2 A^l-^ ^ yswyjF ^3/2 

x,y x,y 



= a^ct ^ P+Uo{T-a, x)-^S{T-a, x) 



-S. . (B.12) 
B.1.1 2— point functions 

To derive explicit expressions for /i, fx and gx in terms of the propagators defined 
above, we use A,B, . . . for the Dirac indices and i,j, ... for the indices in SU{2) 
isospin space. Inserting equations I5.14j) into yields 

/i = -^(^"^" 

x,y,x',y' 

We now perform the Wick contractions of the fields and since the propagator is 
isospin diagonal, each contraction requires the corresponding isospin indices to 
be the same. Due to their Grassmann nature, the sign changes when the order 
of fields is changed 

C'C' CC 5 -5 It^t" 5 -5 = ifTr r'')^ = 

CC' C'C -S- 5 ■ -It^ t" = -iTr fr"")^ = -- 

The Dirac indices give the order of the propagators and Dirac-matrices and we 



E ((C')f(y')75^^|^(C')f (x)C^(y)75^^^Cn,C^(x)) . 



are left with a trace in Dirac and color space 



/i = ^ E (Tr{[C(x)C'(y')]F75[C'(x')C(y)]p75} 



]Bl2l a 



6 



x,y' 

^ ^(Tr{75:s;7575:ST75})^ = ^(Tr^^^^S^)^ . (B.13) 

In the same manner we evaluate the correlation function /p ()5.16|) 

a' 



a 



X 

9 



3L3 

x,u,v 
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Correlation functions 



This is of ttie same structure as for /i and lience tlie only non-trivial Wick 
contraction is [ipC]y [Ci'] with an isospin factor of —3/2 



/p(^o) = ^ E { [^(^)C(u)]p75 [C(v)^(x)]^75} 



x,u,v 



u 



75 



x.v 



u 



2 

X 

^3 



5^(Tr{5(x)5(x)t}) (B.14) 



The correlation function fA{xo) ()5.15j) differs from /p(xo) only by an additional 
7o following '(/'(x). We therefore have 

/a(xo) = ^E(Tr{^(^)^5 [C(v)^(x)]p7o75} 

x,v 



^"^757075 



X 

,3 



These results can also be used to evaluate gp{xo) (j5.18|l if one makes the replace- 
ment 

[^(x), ^(x), C(u), C(v)] ^ [C'(u), C'(v), ^(x), ^(x)] 
in eq. ()B.14|) . The result is 



gpi^o) = ^ E (Tr{ [C'(v)^(x)],75 [^(x)C'(u)]^75 



x,u,v 

.6 



U 



Eini a 



2L3/2 

x,v 

3 



E(Tr{75i?(x)H575i?(a;)7£ 



Again the corresponding axial correlator is obtained by inserting 79. Due to the 
different sign in ()5.17j) we need to put a 70 between ip and the following 75 to 



B.l Summed two-point correlators 
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obtain 

9a{xo) = [^'(^)^(^)]f7o75:R(x)75} 

x,v 



IRToI a- 



^5^(Tr{i?(a;)i?^(x)7o})^. (B.17) 



B.l. 2 3— point functions 

Inserting equations I5.14|l into (j8.15j) yields 

6L6 



X 

,15 



a ^ z6«^^((C')f(u)75^^K.(C')f(v) 



6L6 

x,u,v,w,z 



The only contractions with non- vanishing isospin-factors are those where no fields 
on the same time-slice are paired, namely 

V'C',C^,C'C ^ -SilSknSmj - -|e«%;;.r^,7-„ = -|e"^^TrrV^r'' = |2, 
CC',CV,V'C ^ -5.„,4j^mz - -|e«^%"^.r;',r,; = -|e"^^TrrW = -f^. 

Again the Dirac indices give the order of the propagators and the Dirac-matrices 

-Tr] [C(w)C'(v)]p75 [C'(u)^(x)]^7o [^(x)C(z)]^ 75 



/v(^o)=^E( [^(x)C'(u)],75 [C'(v)C(w)]p75 [C(z)^(x)],7o 



x,u,v 
w,z 



!7 



where we have renamed the arguments of the boundary fields (u ^ v and w ^ 
z) in the second trace under the sum. Using equations ()B.8|) and ()B.10|) the 
propagators from the bulk to the upper boundary can be written in terms of R 
and r\ 

^vM=^$^(Tr{ -i?(x)75 [C'(v)C(w)]f75 [C(z)^(x)] , 70 

+ [C(w)C'(v)]p 7575^(^)^570 [V^(x)C(z)]p75|)^^ . 



x,v 
w,z 
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Correlation functions 



In the next step we insert the expressions for the boundary-to-boundary propa- 
gators St and 5"^ from equations ()B.11|) and ()B.12|) . respectively, and then use 
7575 = 1 in both traces 



4L3 



(Tr {i?(x)755T75 [C(z)^(^ 



u 



Again we use 7575 = 1 and 7570 = —7075 as well as the cyclic property of the 
trace. After inserting the propagators ()B.7j) and ()B.9|) we thus end up with 



/v(a;o) 



4L3/2 



^ (Tr <^ i?(x)75S'r5'(x)^757o + S^'y^Rixy'yo-f^Six 



We now see that the expression is of the form Tr [X + X"^) = 2Re Tr X and the 
final result is therefore 

3 

/v(^o) = ^ (^^Tr {s{xy^,^oRi^HST})^ . (B.18) 



B.1.3 4— point functions 

Preliminaries 

We start by inserting the expression ()4.25p for the improved axial current into 
^.2611 to arrive at 



A^,{x) + acj.dlP'^ix, 



X 



6L6 



Y e'^^^e'^'^^io"^ \Al{x)A\{y) + acKdlP''{x)A\{y) 



^Al{x)acJ)lP\y) + acJ^lP\x)acJ)lP\y) 
= /aa + acAld^JpA + 4'/ap] + a^cl BfdlU^ . 
Here we have introduced the notation 



O 
(B.19) 



fxY{xo,yo)- 



J2e^^'e"^'(p"^X''{x)Y\y)0''), X,Ye{Ao,P}. (B.20) 



x,y 



B.l Summed two-point correlators 
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time 
T 





C' 1=1= C •=C o=C ❖= current /density insertion 




V 

A 




p — p- 




n n 
IL 



a) b) c) d) e) f) g) h) 



Figure B.l: Wick contractions of the correlation function fxYixo^yo)- The gray 
diamonds indicate the insertions of Y and X at times yo £ind xq- 



In order to evaluate these correlation functions in a closed form we define the 
Dirac matrices 

7o75 , if X = y4o [ 7o75 , if ^ = ^0 

, T = < , so that 

75 , if X = P [75 , if Y = P 

X^ix) = tp{x)ElT''ij{x) and = tp{y)TW''ij{y) . (B.21) 



Wick contractions for fxy 

As before the first step is the insertion of the explicit expressions for the boundary 
fields (jOT^IOI|l and the currents/densities (fRTT^ into (jR2n|l . We suppress the 
time arguments of the correlation function to arrive at 

„18 , 

Since we have two currents/densities in the bulk we need to keep track of the 
space-time locations, which we now write as indices. The number of possible 
Wick contractions is significantly larger than for the two- or three-point func- 
tions. 



a) 






CvCw 


/) 


Czc:. 


i^yi^x, 


Ci^y, 


V^xCw) 


b) 


i^xC'u^ C^x, 


Cz^?;, 




9) 




Cz^x, 


Cv^J/, 


^J/Cw; 


c) 








h) 


^.Cu 


a^x, 


Cz^y, 


V^xCw) 


d) 






CvCw 




CzC; 


Ipylpx, 


Ipxi'y, 


CvCw; 


e) 


CzCu' Cv'4'xy 


i'xi'y, 


V^3/Cw 
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Correlation functions 



Isospin factors 

The calculation of the isospin factors corresponding to the contractions a) to i) is 
lengthy but straightforward. With the notation Tr (r" . . . r*^) = Tt a...e the result 
is 



a) 


— 0-^ = 


I n liT' fn<^ n n rt 

~T6 'ij'tl'nt'jn 


1 /7n/^ /'•/if 


Tr 


deba 




= -3/2 


b) 


+5 0^0 = 


1 /in/'' /''fit^ ei fi rt 

16 'ij'jl'nt'tn 


"16^ 


Tr 


daTr 


fee 


= 3/2 


c) 




i 1 ^ahc ^cde —d —CL _6 _e 
^ 16^ 'tj' nV jt'ln 


1 ^dbc^cde 
~ 16^ " 


Tr 


dbTr 


ae 


= -3/2 


d) 




i 1 ^abc ^cde —d —a —b —e 
^ 16 'tj'nl'lt'jn 


1 ^abc^cde 


Tr 


deab 




= 3/2 


e) 


-^in^kj^sl^mt ^ 


i 1 ^abc,cde^d —a—b^e 
16 'nj'jl'lt'tn 


1 ,abc,cde 
~ 16 " " 


Tr 


dabe 




= -3/2 


f) 


_^in^kt^sj^ml ^ 


i 1 ^abc^cde^d ^a^b _e 
^ 16 'nj'tl'jt'ln 


1 ^abc,cde 


Tr 


dbae 




= 3/2 


9) 


_^il^kn^sj^mt ^ 


i 1 ^abc ^cde —d —a _fe _e 
^ 16^ 'ij'nl'jt'tn 


1 ^abc,cde 
~ 16 " 


Tr 


dbea 







h) 




i 1 ^abc ^cde —d —a —b _e 
16^ 'tj'jl'nt'ln 


1 ^abc^cde 
~ 16 " 


Tr 


daeb 












i 1 ,abc,cde—d _a_6_e 
^ 16 'nj'trit'jn 


1 ,abc,cde 
~ 16^ " 


Tr 


deTr 


afe 






Combining this result with the contractions gives the propagators that need to be 
evaluated. Keeping track of the Dirac indices gives the order of the propagators 
in the non-vanishing contractions 

-|Tr {[^..C]p75 [CCw]p75 [Cz^,]pT 
|Tr {[V'.ClpTs [C^.]pS}Tr { [CzV'.lp T [V'.Cwjp 75} 
|Tr { [^,C] p 75 [C^.] p T } Tr { [Cz^.] p H [^.Cw] p 75 } 
fTr {[^,C]f75 [CCw]p75 [Cz^.]pH [^.i^y\^T] 
-|Tr {[CzC]p75 [C^.]pH [V^,.V^,]pT KCw]p75} 
|Tr { [CzC] p 75 [C'Jy] p T K^.] p H [^.Cw] p 75 } . 



a) 

d) 
e) 
f) 



B.l Summed two-point correlators 
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The general correlator fxy 

If we keep the arguments of the fields as indices, the correlation function fxv can 
now be written as 



-Tr {[^.C]p75 [C^.]pH}Tr {[Cz^,]^T [v^,Cw]p75} 
+Tr { [^|JyC] p 75 [C'Jy] ^ T } Tr { [Cz^.] p H [V^.Cw] p 75 } 
-Tr { [ipyC] F 75 [CvCw] p 75 [Czipx] p S [ipxi^y] p T } 
+Tr {[CzC]p75 [C^.]pH [V^.^,]pT KCw]p75} 
-Tr {[CzC]p75 [C;^,]fT [V^.^.jpH [^.Cw]p75} )^ ■ 



Using eqs. (fAl2|l . (|R7jl . (|RTn|l . (|rTi]) and (jRT2|l as well as the properties of 
the Dirac matrices, one arrives at 



fxY = ^ E ( - {R{x)j,STS{y)h5rS{y; x)E} 

"•^ -L3/2Tr {R{x)Rix)h,E} Tr {Siy)h,TSiy)} 
+L'/^Tt {Riy)Riy)^,T}Ti {Sixy^.ESix)} 
+Tr {Riyh,STSixY75ESix;y)T} 

-Tr {5^75:R(x)^75H 75^(1/; x)t75T5(2/)} 

+Tr |5'5.75i?(?/)^75T75S'(a;;?/)^75SS'(x)| . 



If we now specialize to the four cases /pp, /aa, /ap and /pa, the connected 
contributions will combine to form pairs of Hermitian conjugates. 



/pp ^ S = T = 75 



18 



-L3/2Tr {R{x)R{xy}Tr {Siy)S{yy} 
+L^/^Tt {R{y)R{yy}TT {S{x)S{xy} 
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Correlation functions 



AA 



T = 7o75 



AA 



4L3/2 



E 



f 2ReTr {R{x)-f5STS{y)'<-foS{y; x)-fol5} 
-L3/2Tr {Rix)Rix)ho}Tr {S{y)S{y)ho} 
+L3/2Tr {Riy)Riy)ho}Tr {S ix)S (x)^ ^0} 
-2ReTr {R{y)-f5STS{x)^-fo S{x;y)-fol5} 



Zap =^ H = 7o75, T = 75 



AP 



^3/2 



4L3/2 



-2ReTr {R{x)-f,STS{yy S{y; x)jol5} 
+L3/2Tr {;R(a;);R(x)t7o}Tr 
-L'/'Ti {R{y)Riyy}TT {Six)Sixy^o} 
-2ReTr {:R(y)755T5(x)t7o5(x; 2/)75} 



/pa ^ S = 75, T = 7o75 



/] 



PA 



4L3/2 



^B.23) 



;B.24) 



f 2ReTr {R{x)^^STS{y)''^^S{y-x)^^) 
+L3/2Tr {;R(a;);R(x)t}Tr {5'(|/)5(?/)t7o} 
-L^/^Tr {i?(y):R(y)t7o}Tr {5(x)5(x)t} 

+2ReTr {:R(|/)75:Sr:S(x)t^(x; i/)7o75} )^ • (B.25) 

Evaluation of /p^ 

Inserting the expression ()4.25|) for tlie improved axial current into ()8.33|) gives 

7 w+* 



/pa (2/0 +^,2/0) 



6L6 



'^S^ a ^ w(xo) /pA(a;o, Z/o) + acA9^/pp(xo, yo) 
2:0=2/0 



(B.26) 



B.2 Simplifying the correlation functions 
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B.2 Simplifying the correlation functions 

We now specify an explicit (chiral) representation of the Euclidean Dirac matrices 
/ 

7^ = , with the 2x2 matrices Cq = —1, = —icTk , 

where k = 1,2,3 and ak are the Pauli matrices. With 75 = 70717273 and 
P± = |(1 ± 7o) we then have 

If we write the four components explicitly we see that for all matrices X, Y in 
Dirac space we can conclude that 

= XP+ =^ = 2{XP+)a* = {Xai — Xa3, Xa2—Xa4:, Xas—Xai, Xaa — Xa2) 
=^ = Xai—Xa3 = Xa2—Xa4 , (B.27) 

= XP_ ^ = 2{XP^)a* = {Xai+Xa3, Xa2+Xaa, Xai+Xaz-, XA2+XAi) 
=^ = Xai+Xa3 = Xa2+Xa4 , (B.28) 



— P+Y =^ — 2(P+F)*B — {YiB — Y^B, ^23—^4:3, YsB—yiB^ ^43 — ^1 



2B) 



= Y^B-YsB = Y2B-Y^B , (B.29) 

= P_Y = 2(P_F),B = {YiB + Y^B, Y,B + Y,B, Y,b + Y,b, Y^b + Y^bY 
^ = YiB + Y^B = Y2B + Y4B • (B.30) 



From equations ()B.1|) - ()B.5|) we can now derive relations between the Dirac 
components of the summed correlators (only Dirac indices are written explicitly) 



oW^P- ^ Sai + Saz = Sa2 + Saa, (B.31) 



Q^PJS' ^ Q^{s\bH'S%b = {s\bH's\b. (B.32) 
' = ^ RP+ =^ ' ='^' Rai — Ra3 = Ra2 — Raa , (B.33) 
O^P^lt 0^ (r\b-(R\b = (R%b-(R\b, (B.34) 



0^ StP- {ST)Al + iST)A3 = iST)A2+iST)Ai, (B.35) 



O^P^St => {St)ib + {St)3B = {St)2B + {St)ab, (B.36) 



= P.St ^ 75P-5t = ^ P+(75^t) = 

=^ (755*^)15 -(75'S't)3B = (75'S'r)2B-(75'S'T)4i? • (B.37) 
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B.2.1 Explicit form of the correlation functions 



We use A, B, . . . for the Dirac indices and a, /3, . . . for color indices. 

4 3 ^43 



Tr StSI^ 



E E (^r)t {si 

A,B=l a,l3=l 

E E H^t) 

A,B=1 a,/3=l 



T . 



'AB 



BA 

2 Em 



A,B=l a,l3=l 



A=l B=l a, 13=1 



4 E E \(St] 

A,B=1 a, 13=1 



al3 
AB 



:B.381 



The final form for the correlation function /i from ()B.13|1 is therefore 



AB 



\A,B=1 a,/3=l 



(B.39) 



For the more complicated correlation functions it will be useful to introduce a 
modified sum convention, where all Dirac indices that are restricted to {1,2} are 
written in square brackets. In this notation we have 



^B.40) 



The correlation function /v(xo) from ()B.18|) contains the trace of '^^^qR'-^^St- 
Using equation ()B.33|) and ()B.37|) one can see that the Dirac index contracting 
R and 756'^ can be restricted to {1,2} if a factor of 2 is included. Since in this 
representation 75 is equal to the unit matrix if the indices are restricted to {1,2} 
the 75 in front of St can be dropped, i.e. 

■ ■ ■Rab{'1^St)bc ■ ■ ■ = ■ ■ ■'^Ra[b]{St)[b]c ■ ■ ■ 

Similarly, using equations ()B.32|) and ()B.35|) the contracting Dirac index between 
St and S can be restricted to {1,2} if another factor of 2 is included 



{St)ab{s'^)bc ■ 



'^{.St)a[b]{S )[b]c • • 



The result for /v(xo) is then 
2a^ 



L3/2 



5^ Re {S{x)%^i^,^o)BcRixrciD]iST: 



[Dm 



B.2 Simplifying the correlation functions 
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Since the combination S ■ will appear in all fxY as well as fv, we define 



'S{x) {Sry 

The final result for /v(a;o) is then 
2a^ 



Pi 



BID] 



07 

[Am 



L3/2 



^Re I {ST)jD][A]iS{xy)[A]Bh5lo)BcR{x, 



,/37 
'cfDl 



X 



f/ 



/a, /p, 5'a, S'p 



In addition to the above simplifications we also have 



(|R3T|1 A (Ib:32|) • 
(|R33l) A jESil) ^ ...:Rab(^%c 



• • • 'S'ab(5'^)_BC • • • — ■ ■ ■'2SA[B]iS^)[B]C ■ 



■ '^Ra[B]{R^)[B]C ■ 



With these the correlators /p, /a, (^p and (^a can be written as 



/p(a;o) 
/A(a;o) 



13a 
[B]A/ U 
t3a 



B]C I U 



B]A/ U 

X 

^ «^E(fi(-):s^i('5(-)*)l(-'» 



CA 

B]C I U 



(B.41) 



[BA2) 



(B.43) 
(B.44) 



(B.45) 
(B.46) 
(B.47) 
(B.48) 



fxY{xo,yo) 



The expressions for /xy can be written in the same form as fy with the definitions 

mx,yorc[A] = a'j^isi^^yh^mYc^A]^ (^-49) 

(B.50) 



C[A] 



Using these, the expression ()B.22jl for /pp becomes 
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MAU a' 



+2Re {(^t)^^^!^] {S{x;yh,R{y)y;^^^} 
-{x < — ^ y}^^ 



+2Re 



-{x < — > y}^^ 

X 

-{x < — > y}^^ . 

The disconnected parts are built from the propagator products, which also appear 
in /p and gp. However, this is done on each gauge configuration and we therefore 
have to define new correlation functions 



x,y 

X {^y)];^4S{y)XmAio)FD}\ 
[gfWixo, yo) = a'J2{{Ri^)fmWx)%^ciio)cA}^^^^ 

x,y 



x,y 



Without any additional effort these correlation functions can be constructed from 
the estimates for fx and gx in the analysis program. Inserting them into the last 
expression for /pp gives 



B.2 Simplifying the correlation functions 
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/pp = |i^(-Re$:{(p(i/)2f^])*iV5(y,xo)- J 

y 

+Re5^{(p(x)5^JiV5(x,yo)^^, 



u 



-[gf]pp(,xo,yo) + [gf]pp{yo,xo) . (B.51) 



In the expression ()B.23|1 for /aa again two indices in each line can be restricted, 
thus giving an overall factor of 4 



,6 



{R{y)fiB]Wy)%^ciio)cA} 

-{x < — > y} 



L3/2 ^ 

+^'/'{^(2/)^?B](^(l/)^)fB]c(7o)cA 



-{x < — > y} 



u 



y 

-Re ^ I (P{x)l1^^ ilo)cDNo5{x, t/q)?)"^] 



-[5'/]AA(xo,yo) - [gf]AA{yo,xo) . (B.52) 



The evaluation of /ap and /pa from equations ()B.24|) and (|B.25|) is now entirely 
straightforward. However, the antisymmetry under the exchange x ^ y is ob- 
tained only when the two are combined. 
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Correlation functions 



-2Re {(St)^a][b] (^(^)^)Sc ^^^^cd {S{x; 2/)75^(2/))L°^] } )^ 

-2Re { (P(x)g^J (7o)cD (^(x; y)75:R(l/));[^] } )^ 
y 

-Re 5^ { (P(a:)-Jf^,)* (7o)cDiV5(x, 2/o)2,[^]} )^ 

X 

+ [gf]Av{.XQ, yo) + [fi-ZlpAlyo, xo) . (B.53) 



{^(i/)A?B](^(i/)')fB]c(7o)cA} {:^(x)-i,(:s(x)t)£j^} 

+2Re |(^T)fA^[B] (-^la;; Z/)7o75^(z/)) J"^]} )^ 

{Rix)f[B](Ri^)%]A} {^(?/)I>[i.](^(2/)^)g]F(70)FD} 
{^(y)A?B](^(y)')fB]c(70)cA} {S{X)];[E](S{X)%^} 

+2Re { (P(a;)5^ J yhol^mYc^A^} )^ 
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+ Re 5^ { (P(y)g^J (7o)cDiV5(|/,xo)5^]} 
y 



-[5'/]pA(a;o, yo) - [gf]Ap{yo, xo) . (B.54) 
We can now plug the expressions ()B.51|) to ()B.54|) back into ()B.19|) to obtain 

/AA(a;o, yo) = Jaa^xq, yo) + acA[do fpAixo, yo) + 9^/ap(xo, yo)] 
Wc\d^d^fpp{xo,yo) 



2a 



-acAdl 
2a^ 



y 

2a^ 



13/2 \ (^(^)c[A]) (7o)cD No5iy,xo) + acAdoN^{y,xo) 



- 70 

D[A] \ / U 



Nobiy, xo) + acAd^N^iy, xq) 



13/2 1 (^(^)c[A] j (7o)cD No^ix, yo) + acAd^N^ix, yo) 

^ 2a^ 



+acAd^ 



[X 



■ya 



No5{x, yo) + acAd^N^ix, yo) 



C[A] 
•ya 

J D[A] J / U 
•ya 

C[A] I / U 



13/2 ^ 1 ^■^IC\A\ J 

+ [5'/]AA(a;o, yo) - [9f]AAiyo, xo) + acAdo - [fi'/]pA(a;o, yo) - [fi'/]Ap(yo, xo 



+acAdli 



[gf]Ap{xo,yo) + [gf]pAiyo,xo) 



+a^c\dQdl [gf]pp{yo,xo) - [gf]pp{xo,yo) 



(B.55) 



The expressions in square brackets, which multiply the P propagator, can be 
obtained with one additional inversion for xq and yo and each combination of 
external indices [A\ and a. More details about the necessary inversions are given 
in Sections IB.2.21 and IB.2.31 If it is necessary to keep f\p^ as an explicit function 
of ca more inversions are needed. 

The correlation function /pa(z/o + i^yo) (|B.26|1 multiplying the mass term in 
the integrated Ward identity can also be expanded in the same way. 

yo+t 

/pa(2/o + t,yo) = w{xo) /pA(a;o, yo) + CAa9^/pp(xo, yo) 
2a' . J^o+i 



L3/2 



ReY,w{xo)Y,{ [Pi^ 



C[A]) 



.70 



No5{x,yo) 



-ya 



C[A] ) / U 
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+acA 
—cacl 



2a^ 
L372 



yo+t 



xo=yo 



[X 



C[A]J 



C[A] \ / U 



C[A]) 



(7i 



0)CD 



yo+t 



^ w{xo)N5{y,Xo) 



7a 

D[A] 



.xo=yo 



2a^ 
I372 



)7a 
C[A] 



A^5(z/,a;o)c"A] 



xo=yo 



yo+t 



a;o=j/o 



^o( - [^/]pp(a;o,yo) + [gf]pp{yo,xo] 



yo+t 



a^w{xo)( - [gf]pA{xo,yo) - [9f]Ap{yo, 



(B.56) 



xo=yo 



B.2.2 Sources for the inversion of the Dirac operator 

With the definition of the propagator S{x; y) as the inverse of the Dirac operator^ 

y)^^ ■ S{y; zf^^ = a~'^6^Ja^6Ac , (B.57) 
one can easily calculate the action of D on the summed correlators. 

Em cta^ 



D{x,yrA%S{yr 



B[C] 



L3/2 



^ D{x, yTlsSiy. a, z)f^[f/o(0, 7.r%,[P^]n 



[c] 



L3/2 



Ct 



5a,xo [f/o(0,x)"l]^^ [p_ 



Ja[c] , 



D{x,y)f^R{yr< 



Ct 



)[C] 



:B.581 



(B.59) 



The right-hand sides of equations ()B.58|) and ()B.59|) are therefore the sources 
to be used in the calculation of S{x) and R{x). For all index combinations of 
these propagators that were used, six inversions, corresponding to the possible 
combinations of 7 and [C], are needed. Similarly we obtain for the expressions 
containing the full propagator S{x] y) 



C[A] 



5a 
E[A] 



a ^5^0X0^ BE [R{z)] 



f3a 

E[A] ■ 



(B.60) 



^Here D stands symbolically for D + 6D+rno from (|A.l 
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Thus, to calculate and A^os, the knowledge of the full propagator S{x; y) is not 
necessary. Instead, one uses R as the source for an inversion to obtain the desired 
product of propagators. When is used with a time derive, this is included in 
the source, such that e.g. d^N^ly, xq) also requires only one inversion. The same 
applies to the correlator 

yo+t 

xo=yo 

appearing in both the bare and the improvement term of /p^. 



B.2. 3 Counting inversions 

Here we give the structure of the program that calculates the correlation function 
necessary for the estimate of Za- We write /^^ as f AAl+CAf AA2+c^f AA3 and /p^ 
as f PAl+CAf PA2. On a given gauge configuration the program proceeds according 
to 



for all 7 and [C] <^ Solve for and 

accumulate result into /p, /a, St, fi and P. 



for all 7 and [C] < Solve for i?(T<r)*pj and 

accumulate result into gp, gx, fy- 

Solve for NQ^{-k,Xo)*Jij^ and 
accumulate result into f AAl and f AA2. 

Solve for c?q A^5(*, xo)*!!^.] and 
accumulate result into f AA2 and f AA3. 

Solve for iVo5(*, Z/o)tj!c] 

accumulate result into f AAl, f AA2 and fPAl. 

Solve for (9q A^5(T»r, yo)tjlc] ^^d 
accumulate result into f AA2, f AA3 and f PA2. 

Solve for Efo=.o^(^' ^o)?^ 
accumulate result into f PAl and fPA2. 

The "accumulate into..." corresponds to the trace over color and (restricted) 
Dirac indices in the correlation functions. For the connected correlators this is 
done by taking the scalar product of P with the result of the inversion, taking 
into account temporal derivatives or sums and the correct combination of Dirac 
indices as given by ()B.55jl and ()B.56jl . 
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Thus, a calculation of /a, /p and /i needs six inversions on each gauge con- 
figuration, while six more are required for gp^, g-p and /y. A total of 36 is needed 
for the correlation function f\p^ and only six more for the volume correlator /p^. 
However, even this is negligible compared to the cost of the inversions in the 
molecular dynamics. One should note that if the value of ca were fixed at run- 
time, two of the additional inversions could be saved by using the linearity of the 
Dirac equation, e.g. by solving 



D{z, vYbc No5{y, xo) + acAd^N5{y, xq) 



■ya 



C[A] 



-^Szo,xo-a{l5)BE [R{z 



E[A] 



thus reducing the total number of inversions from 42 to 30. An overview is given 
in table Table Ell 



correlation functions 


^ inversions 


/a, If, fi 

/a, /p, /i, fl'A, gp, /v 


6 

12 


/a, /p, /i, gA, gp, /v, /aa 

/a, fp, fi, gA, gp, /v, Jaa, /pa 


36 (24) 
42 (30) 



Table B.l: Number of inversions needed to compute the correlation functions. The 
numbers in parenthesis refer to the case, where the value of ca is fixed at run-time. 



Appendix C 

Transforming the integrated 
Ward identity 



We start by isolating the contact term in the volume integration in ()8.29j) 

r'"yo+t 
ya-t 



dV yd^xy ^ dxoe'^''^(P'^(xo,x)A|;(2/o,y)a: 

/r ryo+t , 

dV /d^x / dxoe"'''=(P'^(xo,x)A|;(yo,y)a 

dV / d^x / dxo e"'7 P'^lxo, x)A^(yo, y)axt 
J Jyo-<^ 

abc I T->a 



-2m Jd'y Jd'^J ^ dxoe'^^^^P'^(xo,x)A^(yo,y)axt; • (C.l) 

We will now use the partial conservation of the axial current to relate the two 
contributions from the surface integral in eq. ()8.29|) . Using Jd^x9fc/fc(x) = 
(for periodic spatial boundary conditions) the partial conservation reads as an 
operator identity 



'dV doAliy) = J dV d.Aliy) = 2m J d'y P\y) 

/r ryo-<^ r 

dV^^(2/o-t,y) = /dV4(|/o-e,y)-2m / dxo d^yP\xo,y) . (C.2) 



'yo-t 

We will also use ()(I2|1 in the form 

/r ryo+t p 

d=^xAS(yo,x) = j d3xAS(yo+t,x) - 2my dxo j d=^xP"(xo,x) . (C.3) 

These relations can be used only in matrix elements with fields that are not 
located in the integration region since otherwise the axial variation of the latter 
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Transforming the integrated Ward identity 



will appear as an additional term. We therefore have to be careful how far we 
can shift the current insertions. Due to the antisymmetry of e"'"^ the second 
expression in the first line of equation ()8.29j) can be rewritten as 



Jd'y y"d3xe'^^'^(A^(|/o,x)A^(i/o-t,y)axt 



dV /"d3xe'^*'=(AS(yo,x)A^(yo-e,y)a 



EH 



\--\j\c<J7--/--\j\^'^ - -I ■u / - CXt 

-2m /dV /d'x dxoe'^^^(ylS(yo,x)P^(xo,y)axt 



dV yd3xe'^^'=(Ag(2/o+t,x)A^(yo-e,y)a: 

dV yd=^x y dxoe'^''^(P'^(xo,x)A[;(|/o-e,y)a 
+2m ydV y"d=^x y"'°^dxoe'^''^(p'^(xo,y)A|;(yo,x)axt) • (C.4) 



The last term cancels the last integral in the split up integration (jCIlj) and in 
the limit e — > the first combines with the result from the integration over the 
upper surface. Thus, inserting ()(].4jl and (jC^.ljl into equation ()8.29j) yields 



lini 2 /d'y /d'xf'"-(AS(!,„+t,x)/lS(!,„,y)a 

dV yd^x y dxoe'^^^(P'^(xo,x)A|;(yo-e,y)a 
dV / d=^x / dxo e»''^(p»(xo, x)A|;(yo, y)axt 
-2m y"dV yd^x y"'°'"'dxoe'^^'=(p'^(xo,x)A|;(yo,y)axt) .(C.5) 



In the third line we can replace A\^{yQ—e) by Ao(?/o) if "we start the time integration 
at Uo + e. This preserves the order of the operator insertions and allows us to 
combine two of the time integrations with the result given by 
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Jim Jd'y y"d3xe"^^(A^(|/o+t,x)4(yo,y)a 

-m Jd^y Jd^^ y"'°^'dxoe«''^(p'^(xo,x)A^(yo,y)axt) . (C.6) 



The contact terms are the same as those appearing in ()8.28j) and thus integrable 
under a four-dimensional integration. The hmit of e — > can hence be performed 
and the final result is 



2 Jd'y (V,%yo,y)0,: 

fd'y jd^^e-'^(^Al{y,+t,^)Al{y,,y)0, 

/r ryo+t , 
dV yd=^x J dxoe'^''^(P"(xo,x)4(yo,y)a 



Appendix D 

List of simulation parameters and 
results 

In Table ID. II the simulation results for and are collected. The algorithm 
is specified in the same manner as in Table 16.11 The number of measurements 
^meas IS multiplied with the number of replica and as before Tmcas is the molec- 
ular dynamics time between consecutive measurements. The PCAC mass am is 
defined through ()8.37j) and the column Zy refers to the definition ()8.3(i|l before 
taking the chiral limit. 

For completeness the results for Za are also given for the estimate using 
only connected diagrams in the evaluation of the correlation functions /xy? see 
Section [8.3.41 This is denoted by Z'^'^ . An (additional) superscript "old" refers 
to the estimate from the massless condition from ref. [11], where the vnf\p^ term 
is neglected. 
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algo (3 K L T Nmeas Tmeas 


am Za ^a°° ^a^'^ ^™n,oid 


Ha 5.200 0.13550 8 18 16-200 4 
Pi4o 5.200 0.13550 8 18 16-40 10 
Pi4o 5.200 0.13560 8 18 16-225 3 
Pi4o 5.200 0.13570 8 18 16-230 2 
Pi4o 5.200 0.13580 8 18 16-230 2 


0.01718(90) 0.7301(173) 0.8411(80) 0.5039(159) 1.0124(224) 0.7509(6) 
0.0159(11) 0.7186(295) 0.8455(108) 0.5026(156) 1.0174(227) 0.7497(14) 
0.01310(68) 0.7157(137) 0.8212(96) 0.5546(117) 0.9627(123) 0.7471(7) 
0.0088(11) 0.7134(126) 0.8302(70) 0.6222(149) 0.9114(113) 0.7447(8) 
0.00194(81) 0.7176(114) 0.8588(99) 0.7032(115) 0.8773(104) 0.7424(14) 


Pi4o 5.290 0.13625 8 18 16-50 2 
Pi40 5.290 0.13641 8 1 8 1 6-120 2 


0.0031(18) 0.7527(102) 0.8103(167) 0.7391(108) 0.8437(179) 0.7507(19) 
-0.00512(61) 0.7540(124) 0.8378(73) 0.7421(120) 0.8081(62) 0.7490(12) 


Ha 5.500 0.13606 12 27 16-25 6 
Ha 5.500 0.13650 12 27 16-44 3 
Ha 5.500 0.13672 12 27 16-80 3 
Ha 5.500 0.13672 8 18 1-318 4 


0.02254(26) 0.8417(222) 0.8077(26) 0.3918(137) 1.0732(261) 0.7853(14) 
0.00758(27) 0.7987(153) 0.8100(45) 0.6063(143) 0.8630(57) 0.7738(8) 
0.00041(25) 0.7888(32) 0.8048(54) 0.7861(33) 0.8096(56) 0.7650(21) 
-0.00168(62) 0.8105(64) 0.8168(38) 0.8148(58) 0.8111(42) 0.7750(45) 


Ha 5.715 0.13665 16 36 1-106 2 
Ha 5.715 0.13670 16 36 1-54 2 
Ha 5.715 0.13670 12 27 4-62 2 


0.00194(57) 0.8142(135) 0.8079(31) 0.7811(60) 0.8199(20) 0.7827(11) 
-0.00060(69) 0.8004(26) 0.8120(30) 0.8014(23) 0.8098(28) 0.7793(20) 
-0.00100(34) 0.8021(38) 0.8182(18) 0.8071(34) 0.8138(19) 0.7861(18) 


Ha 7.200 0.13420 8 18 1-220 2 
Ha 7.200 0.13424 8 18 1-164 2 
Ha 7.200 0.13424 12 27 16-50 2 
Ha 7.200 0.13424 16 36 1-80 2 


0.00029(45) 0.8721(24) 0.8772(18) 0.8699(22) 0.8787(18) 0.8573(9) 
-0.00028(42) 0.8683(22) 0.8732(18) 0.8707(36) 0.8716(28) 0.8553(6) 
-0.00049(15) 0.8685(23) 0.8717(8) 0.8717(23) 0.8697(7) 0.8543(18) 
-0.00023(41) 0.8678(18) 0.8670(17) 0.8694(14) 0.8652(11) 0.8508(18) 


Ha 8.400 0.13258 8 18 4-40 2 
Ha 8.400 0.13262 8 18 4-45 2 


0.00023(40) 0.8990(28) 0.8956(16) 0.8962(26) 0.8973(19) 0.8839(15) 
-0.00183(42) 0.8998(25) 0.8953(13) 0.9184(39) 0.8845(15) 0.8826(7) 


Ha 9.600 0.13140 8 18 4-100 2 
Ha 9.600 0.13142 8 18 4-125 2 


0.00021(15) 0.9137(14) 0.9154(7) 0.9110(18) 0.9166(8) 0.9040(4) 
-0.00059(15) 0.9118(12) 0.9155(7) 0.9188(18) 0.9121(10) 0.9034(4) 



Table D.l: Summary of simulation parameters and results for and Zy. 
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